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ABSTRACT 


Tke primary intent of this work is to investigate 
computer aided compensator deSign using classical 
freguency respcnse technigues in conjunction with the 
technigues of modern mathematical prcgramming. A 
computer program to perform the automated design task 
is presented. This particular algorithm is based on 
the ccnstrained optimization technique introduced by 
M. J. ECX. 


In particular, the desired cpen lcop frequency 
response is specified for a number of discrete 
frequency feints over the frequency range of interest. 
Then tke ginimization routine is used to vary the 
compensatcr parameters in such a manner as to minimize 
a cost functional based on the difference between the 
actual and desired open loop frequency response of the 
compensated system. To illustrate the algcerithn 
several detailed examples uSing the program are 
Fresented. 
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I. INTRODUCTION 


Concurrent with the proliferation cf the digital 
ccmputer as a viable, indeed often indispensable, 
engineering tool, the development of state-space methods 
resulted in the application of the digital computer in the 
area of modern control theory over the past decade anda 
half. Initially this application was primarily cf an 
analysis nature, but as the theory of optinal control came 
into wider publicaticn and acceptance, the role of the 
digital ccmputer began to take on aspects of the cesign 
functicn. This situation is graphically illustrated by the 
myriad of industrial proprietary and university develcped 
ccmputer programs for computer aided control system design 
based ufcn state variable and optimal ccntrol thecries. 
However, as cointed out by Mitchell [11] comparatively 
little work has been done to extend and apply the 
ccmputaticnal fower, speed, and accuracy of modern digital 
ccmputers tc the proven classical theory and methcds of 
ccntrol system design. After this essential abandonment of 
the classical sethods (particularly by the academic and 
thecretical factions of the control engineering ccmmunity) 
in favor of the modern theory and techniques, recent 
puklicaticns indicate a realization of the need that both 
approaches to system design be used in a complimentary 
rather than sutually exclusive manner. In particular the 
works of Ccffey, MacFarlane, Mitchell, Page, Rosentrock, 
Stear, [5,8,11,14] and others show a renewed interest in the 
tried and prcven classical thecry cf design coupled with the 
new technigues and approaches developed by proponents of 
mcdern ccntrcl theory, as applicable to the develcrpment cf 
usartle ccmputer aided control system design algorithms. 
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Centrcl systen design based on frequency domain 
specifications, for linear time invariant parameter systems 
using classical methods, has traditionally been done 
essentially using trial and error technigues. That is, cnce 
the form and location of the compensator has been decided 
upcn, the design process reduces to cne of iterative 
analysis until the specifications are wet or it beccmes 
obvious tc tke designer that the compensatcr form chosen 
will not satisfy the specifications. When systems cf any 
complexity are invclved this method quite cften transforms 


intc an extremely tedious and time consuming process. 


The availability of a good analysis pregram coufled with 
a digital ccmputer operating in an interactive mode, to 
prcvide sinimal turn around time, can be used to somewhat 
Binimize the tedium of this type of design procedure. By 
varying the ccmpensator parameters the designer may then 
observe the effects upon system performance withovt having 
to expend ccnsiderable effort in repeated calculaticns. 
This design by iterative analysis may be ccntinued allowing 
the cosputer to perform the calculations until the designer 
is satisfied with the respcnse obtained. With the recent 
refinements in mathematical programming techniques, it is 
possible to relegate some of the design decisions directly 
to the ccaputer prcgram, thus even further minimizing the 
tiwe and efforts reguired of the design engineer. The 
ultimate coal, which is of course at present a long way down 
the road, wceuld be to input only the plant dynamics ard the 
desired specifications and have the computer perforgz the 
design process completely, without any interaction with the 
design engineer. 


In this thesis an investigation of the autcmaticn cf the 
design precedure is undertaken. In particular, giver the 
desired cpen loop response in the frequency domain a 
finimization routine is used in conjunction with an analysis 
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algcrithe tc systematically vary the compensator parameters 
tc approximate the desired response in acccrdance with a 
performance measure which indicates the deviation cf the 
actual open locp frequency response from the desired open 


iccp frequency response. 
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A. GENERAL AUTOMATED DESIGN PHILOSOPHY 


The word design, as used in many engineering texts of 
various disciplines, runs the gamut of describing a process 
dealing essentially with scaling already existing 
configurations and components to one dealing with a kighly 
creative activity allowing a maximum amount of freedcm in 
satisfying a set of requirements or specifications. Ona 
Gay to cay Easis the majority of the design work 
acccmplisted ty eéngineers lies somewhere Letween these two 
extremes, ccmbining elements of relatively reutine, straight 
fcrward calculations with those of creativity and 
imaginaticn. 


The cverall scope of the design process in general is 
indeed sc fcrmidable that it often defies complete 
descripticn in other than general terms, much less 
guantificaticn in precise mathematical language. This 
latter is unfortunately a necessity in achievine the 
autcmaticn cf this process within the limitations of the 
presently available hardware and software accessable Ly the 
engineer. Largely for this reason, when considerinc the 
automatic design problem it must by necessity be viewed, for 
the near future at least, in terms of a limited scope. That 
is, such items as creativity, experience, and intuition, 
which are undoubtedly elements of the total design process 
aS viewed ir its entirety, cannot be achieved by a cconzputer 
as the device exists at the present tine. 
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One cf tke key elements of any design oprcecess is 
decision making. In the overall design process the engineer 
makes decisions based upon the information available to hin 
at the time cf the design undertaking. Quite often these 
decisions are based on a comparison of numerical values of 
various elements of the design configuration. It is ir this 
specific area of the overall design process that the rodern 
Gay computer can serve as a Significant tool for the 
encineer invwclved in design, particularly when a large 
number cf ccapariscns (expressable in numerical terms) are 
invclved. In other words, the computer can be used ina 
ccmkinaticn analysis and elementary decisicn paking role to 
aid the engineer in the overall design precedure. Ey means 
of analysis cumerical values may be obtained, forming a 
basis frcem which comparisons can be drawn. The analysis 
portion is generally essential for exact design and forms 
the basis fer any optimization. The speed available in the 
digital ccmputer allows many calculations and comfariscns, 
over a wide range of possible solutions, to be performed 
which would Ee prohibitive if peformed by hand. 


The employment of the digital computer as an aid in the 
aggregate design process is perhaps as much an art ag the 
design process itself. While the use cf computer aided 
design prcgrams can free the engineer from many burdensome 
and time ccnsuming calculations, the individual enogineer 
must still previde the framework in which tc interpret the 
results in the hard light of physical significance. The 
engineer must also be capable and willing to define the 
Frcblem in a form with which the computer can work. More 
often than net, this latter requires a much sharper, clearer 
definiticn (and subsequently understanding) of the problen 
than if the same problem were teing designed long hand. 
Unlike the analysis problem, which generally possesses a 
specific answer, the design problem may not have a eunigue 
scluticn which further complicates any computer aided design 
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algcrithz. Ancther aspect that should not be overlocked is 
the fact that when using a computer aided design program the 
engineer cften receives, if nothing else, an indicaticn of 
the soluticns which are not feasible. This in itself is 
often nct a trivial result especially in dealing with 
systems cf any complexity. 


Ancther important factor which should nct be overlooked 
is that cnce an answer is obtained using a computer aided 
design algcrithm scme type of sensitivity study will more 
than likely ke necessary. This is largely due to the fact 
that ccmpcnent values in phySical hardware cannot generally 
be contrclled to the precision with which a computer does 
its calculaticns. Therefore, while a particular algcrithno 
for use in automated design may generate what is indeed a 
theoretically flausible solution the implementation of such 
a soluticn gay not be possible if it is extremely sensitive 
te the parameter values involved. 


B. SOME ASFECTS CF COMPUTERIZING THE CLASSICAL DESIGN CF 
CCNTROL SYSTEMS 


Classically the purpose of the design problem in 
feedback contrcl systems has been to make a given plant 
satisfy a set cf performance specifications imposed upen the 
system. fLesign of linear, time invariant ccmpensators for 
the ccntrcl cf linear systems has formed a major pcertion of 
the design effcrts of conventional servomechanisms. In the 
process cf investigating the automation of this design 
precedure, certain key elements were found to be generally 
necessary in istpdementing a basic practical ccsputer 
algcriths for autorated design, regardless cf whether it is 
based on time cr frequency domain approaches. These 
elements consisted of: 


18 


1. A ccmglete description of the plant 
cr process to be controlled. 


2. Specifications describing the 
response desired from the system. 


3. Some type of criterion functicn or 

Be Cao tiaeaetteadecrec; en 

4. The form or structure of the 

tr epsacor Or ccntrol to be 
These items appear in varying degrees aS a common 
denominatcr amcng present day automated design algorithms 
{1,5,11,13] regardless of the specific methods empfioyed or 
systems keing considered. That is to say, these elements 
appear in algorithms designed for time demain, trarsforn 
dowain, stochastic frocesses, and nonlinear systems. 
Generally items 1. and 4. are employed in some fcecrm of 
analysis algcrithm while items 2. and 3. form the basis for 
a Minimization (cr maximization) routine that prevides a 
foundaticn fcr some type of decision making withir the 
algorithn. This idea of using a performance measure upon 
which to kasé scme type of decision making routine is a 
direct carrycver from develcpments in the computerized 
design of oftimal control systems which relies heavily upon 
the idea cf minizizing some criterion function. Unlike the 
technigues cf cptimal control however, the cost functicnal 
as used in this ccntext is nct necessarily a function cf the 
states and tke applied control, rut rather a measure of the 
deviation of the actual response from some desired response 
which has been specified by the designer. This is ferhags a 
fine distincticn but one that must be recognized and kept in 
Bind. One of the sany consequences of this is that there is 
no ilicnger a strcng motivation to restrict the performance 
measure tc a guadratic form in order that the second sethod 
of Liapunov may te used to frove that a stable system will 
result. Rather, the form of the cost functicn in 
ccmputerized design based upon classical ccntrol thecry is 
Chosen tc have properties which aid in the convergence of 
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tke numerical pinimization routine used in the algcrithm or 
which has scre inherent physical significance in the design 


prcecess. 


C. A FCEMULIATICN OF THE AUTCMATIC DESIGN ERCBLEM IN THE 
FREQUENCY DCKAIN 


In this thesis the automatic design prcblem based on 
theory develcped for classical design methods using 
freguency dczain procedures is investigated. Several 
approaches Fased on frequency domain techniques have keen 
presented by various people to automate the control system 
design precess [5,11,13,14,20]. These algorithms range fron 
methods cesigned primarily for use with multivariable 
systems to specialized interpretations of frequency domain 
data such as vector frequency domain technigues and inverse 
Nyguist setkcds. The suitability of these various methods 
have keen demcnstrated by their respective authors and 
others scrking in the area of computer aided control system 
design. While each represents a step forward in the line of 
progress tcward automating the design process using 
frequency dcmain techniques the Suitability of each 
algorithm is dictated to some extent by the particular type 
of design prceblem being solved. Many cf these methods 
require gradient calculations, and those that do not rely 
primarily on simplex nonlinear programming algorithms. The 
former methcds restrict the ferm of the cost functicn to be 
Ririmized in the sense that it must be differentiakle and 
prcblems inherent with numerical differentiating algorithms 
are obvicusly involved. The latter methods invclving 
Sinplex algcrithms have difficulty in handling constrained 
tHinimizaticn problems, particularly when solutions exist 
hear or cn the constraint boundaries. 
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One cf tke advantages of the frequency response apfroach 
lies in the fact that it lends itself easily to the use of 
cpen-loor frequency response data in designing closed loop 
systems. The use of frequency response methods provides a 
graphical pertrayal of the system response that lends itself 
well to quick interpretation Ey the designer. Also since 
approaches such as those developed by Bode were intended to 
give Ssitrle approximate methods for use in hand 
calculaticns, adaptation of these methods in a computer 
algoriths allcws the engineer te plan his design using such 
technigues and then employ the aid of the ccmputer without 
the requiregzent cf significant problem restructuring in 
order to cast the problem in some particular format not 
generally familiar to the designer, but necessary for 
solution Ey a computer algorithn. With systems invclving 
ccmplex transfer functions such restructuring can itself 
become quite kurdensome thus defeating one cf the primary 
reascns for turning to computer aided design. Perhaps it 
shculd be peinted cut at this fpoint that even using a 
ccuputer algorithm to aid in the design process does not 
totally relieve the engineer of all hand calculations. Time 
must be taken to carefully map out a plan of attack cna 
particular freblem. In essence the designer must thrcughly 
understand the problem and this is generally achieved by 
scgé preliginary rcugh sketches and calculations before the 


ccmputer ever sakes its appearance in the design process. 


Throughout this thesis, the vehicle used for presenting 
the autcmated design algorithm is the Singleée-input 
Single-output unity feedback system of the form shown in 
figure II-1. 


2 


General System, Block Diagram 


figure II-1 
The plant is considered to be linear, time invariant, and 
ccapletely kncwn. Thus the plant may ke descrited Ly a 
rational transfer function in the Laplace variable s of the 


form _ 
Tle + 2) 
i= 
ae m-(s + p.) 
1=0 
where 46 and P. may be complex constants. The compensator 


G (s) is alse considered to be a rational transfer furction 
c 


an s. 


The specifications are given in terms cf the system open 
locp freguercy response, and the desired respcnse is 
achieved by reshaping the system open loop frequency 
response thrcugh the use of the compensator in the _ feed 
forward fath. This basic procedure is well known tc most 
control engineers and is by no means a novel approach to the 
design roklem [3,.10,17,18], but through the years it has 
found extensive use in the practical design of centrol 
systems. 


ae 


The methcd proposed and under investigation here is a 
relatively straight forward adaptation of this technigue 
with a ccuputer alcgorithm harnessed to do the ccmpersator 
parameter variaticn in order to systematically reshare the 
system open locp frequency response with a pinimnum amount of 
calculaticn and iteration required on the part cf the 
designer. The desired open loop response is specified by 
means of a sagnitude and phase profile cver the frequency 
rarge of interest. This is nothing more than a set of 
desired magnitude and corresponding phase values at discrete 
frequency foints selected by the designer. In addition a 
specific form or structure is postulated for the compensator 


transfer function G (s). At this point cf the autcmated 
c 


design process the algorithm relies heavily cn the 


experience and irtuition of the individual designer in 
checsing a structure or form of compensatcr that will have 
scmeé reascnakle prcbability of success in resShafinc the 
frequency respense curves in order to satisfy the design 
specificaticn. Since no firm analytical kasis exists for 
the deterrination of a Specific structure for the 
ccmpensatcr, the computer algorithm must be given this 
informaticn or an iterative method based on some type of 
criterion must be used in order to estaklish the form. This 
latter wethcd was felt to ke much too time consuming and 
restrictive, because of the wide range of variations in 
Fessible structures that might have to ke considered, and 
essentially reyond the scope of this investigation. Thus 
the choice of the compensator form is relegated to the 
program user as a design decision that must be undertaken 
Prior tc using the algorithm to determine specific values 
for the ccefficients of the compensator transfer function 
that will yield the desired response. This procedure should 
nct present a major stumbling block for the designer 
prcvided that he possesss a thorough understanding cf the 
prceklem tc re solved. In fact, by inputing the specific fcrn 
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of the ccmpensator the designer may investigate compersator 
configuraticrs which might otherwise be considered 
unfeasible fer the solution of the problem if some type of 
iterative algcritha were used in determining the structure. 
With these items ferming the basic inputs tc the algcrithm a 
ncnlinear perfcrmance measure coupled with a minimization 
reutine is tken used to minimize the difference between the 
desired and actual frequency response of the ofen loorp 
transfer function given by, 
T(s) = ee oe 


The perfcrmance measure (or cost function) selected is cf a 
normalized form frimarily as a matter of convenience in 
interpreting the algorithm's results regardless of the 
particular yreblem being considered. That 1s , the cost 


function -tas a general form given by 


A (jx) 
D (jy) 
where D(jw) represents the desired open loop frequency 


J(jw) = £(1.0 - 


response cf the entire system, A(jw) is the actual ofen loop 
response, ard w represents the set of discrete freguency 
values at which the cost function is to be evaluated. As 
mentioned previously this general form cf normalized cost 
function was selected primarily as a matter of convenience 
in that a farticular cost function can be constructed in 
such a manner that the optimal or smallest value cf J (jw) 


under ideal scluticn cf the problem will be zero. 


By inputing the desired response in the form of the 
magnitude and phase prefile described above, ncrmal 
freguency dcnzain specifications such as gain margin, fhase 
Margin, and tandwidth (here defined as the open loorg gain 
crossover frequency) may be supplied to the algorithm in one 
type of fcrmat. With these specifications forming part of 


the gair ard phase profiles they are automatically 
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inccrporated into the overall cost functicn each tise the 
algorithm evaluates the frequency response cf the ofen loop 
System cver the range of frequencies specified ky the 
designer, thus also eliminating the requirement of 
specialized ccmputation within the algorithm to check for 
satisfacticn cf these particular specifications as the 


compensatcr farameters are varied. 
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A. GENEBAL 


A variety of functional minimization techniques have 
ccue into use in recent years are discussed in ref. [8]. 
The minimization algorithm used in this work was developed 
originally ty M. J. Eox, [4], and is capable of finding the 
hinimum cf a general nonlinear cost functicn, composed of 
several variatles, within a ccnstrained regicn. The ccmplex 
method of Eox has reen implemented in program form as _ fart 
of the stancard subroutine library at the Naval Postgraduate 
School W&. EF. Church Computer Center. The subrcutine 
(BCXPLX), was originally programmed by R. R. Hilleary cf the 
staff of tke Naval Postgraduate School Computer Center and 
has subsequently been modified slightly Ey the author in 
order to render it more suitable for use in the scheme of 


the compensatcr optimization progran. 


In additicn to its availability as a standard 
subroutine, several other features of this fparticular 
finimizaticn algorithm made it an attractive choice fcr use 
in minimizing a cost function based on frequency domain 
specifications. While admittedly more efficient and 
sophisticated algorithms for function minimization are 
available, tke flexibility and ease of programming of the 
cctplex seetkcd of Bcx are among considerations that should 
not be overlccked. This method does net require that the 
derivatives cf the cost function or implicit constraints be 
calculated, thus avoiding the sundry prektlems typically 
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asscciated tsith numerical methods for the calculaticn of 
derivatives, and eliminating the restriction that the 
derivative cf the cost function exist over the range of 
interest. In implementing the minimization method the cost 
function need not be expressed explicitly in terms cf the 
variables tc be manipulated by the program in the 
Birimizaticn process. Unlike many of the cther ninimization 
prcgrams avéilable, BOXPLX, also provides a restart 
capability. That is, once a possible minimum value cf the 
cost functicnal is reached the program automatically 
restarts frem rardomly generated initial values cf the 
variables. f&khile this method by no means guarantees tkat a 
glcbal einimum will be achieved for every conceivable 
prcktlem the restart capability does at least sorewhat 
attempt tc avoid the problem cf the minimization terminating 
on a local winima within the search area. This methed is 
also tc scme extent scale independent in that the size of 
the initial cecmetric figure used in the search for the 
Binimum value cf the objective function is rcughly scaled to 
tke crder cf the problem variables, through the use of the 
difference ketween the lower and upper bounds cn the 
variables. 


As is the case with other ginimization algorithms, this 
Particular cre is not void of inherent disadvantages and 
difficulties. Present in any Mminimizaticn algorithm of 
course is the recurring problem that no guarantee can be 
made of achieving a global minimum fer all classes of 
problems. The ccmplex method cannot handle eguality 
constraints without modification cf the algcrithm. Also as 
the number of variables increases the method rapidly becomes 
inefficiert. Finally the unconstrained probkleng is 
apparently mcre efficiently handled via gradient technigues. 


B. DESCKIFIICN OF THE ALGORITHM 
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The siniwization method developed by Ecx has its kEasis 
in the seguential simplex techniques of linear prograaming, 
but is specifically designed to obtain solutions to freblems 
cf ccenstrained optimization and to avoid such prokless as 
the inakility of ccnstructing the first siuaplex. Thus the 
resulting scluticn of the minimization problem is a set of 
variables wkich may lie at the extreme edges or boundaries 
cf their permissible ranges. The problems generally 
associated with ccnstraints on the variables is essentially 
attacked Ey tsing a flexible geometric figure having n+l, or 
more, vertices and capable of expanding or contracting, in 
any or all directions. Therefore, unlike many of the 
Sinplex alccrithms no provision is made to maintain a 
regular cecnetry between the various ~ points in an 
n-dimensicnal space, but rather a changing structure (dukbed 
a ccmplex), capable of flattening, rounding corners, and 
eventually ccllapsing on the point representing the ninimun 
of the oktective function, is employed rendering the sethod 
more flexible in handling constrained minimization froktlems. 
The numker cf independent variables in the cbhjective 
function and ccnstraint equations dictate either directly or 
implicitly tke dimensionality cf the space tc be searched. 
The vertices saking up the complex are continually rejected 
and generatec as a search for the minimum cf the okjective 
functicn is carried out with each new vertex generated being 
reguired to satisfy all the imposed constraints. 


The ccnmplex method searches for the minimum of an 
objective function J(x), where the vector x represents a set 
of n variatles, within a region of space Ecunded Ly ufrer 
and lower linits on the variables (explicit constraints) 
given by, 


boeeex S 0 , 1=1,2,3;,---,0 
1 al ak 
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and in addition subject to the restrictions imposed by the 
ccnstraint furctions (implicit constraints) cf the forn, 
a, @ 2 0 Jz1, 2737000 7Me 


As with tke simplex methods this is an iterative algorithn. 
It uses a set of k 2 n+1 points which form the vertices of 
the complex and simultaneously satisfy all the isposed 
ccnstraints. Initially only one pcint or vertex oust 
satisfy all the comstraints and serves as a means of 
generatirg the other k-1 vertices through the use of 


pseudorandce rumbers r., uniformly distrirkuted cver the 
al 


clesed interval from zero to one. The additional k-1 


vertices cf the first complex are established by the 
fclicwing fcruula: 


ad 
i) 
tH 


wat tea. = Le) 
1 ah i a a 

As can re seen frem the formula, the vertices generated in 
this fashion will always satisfy the explicit constraints, 
however, €ach cne gsust be checked to insure compliance with 
the implicit ccnstraints. If an implicit constraint 
viclation is indeed found to cccur at this stage then the 
trial vertex is simply moved, repeatedly if necessary, 
halfway in tcward the centroid of the cther already accepted 
vertices, until ultimately a feasible eint is fcund. 
Continued application of this precedure will result in the 
gceneraticn cf the k-1 additicnal points required fcr the 
initial ccmplex within the feasible region defined by the 
ccnstraints. 


Cnce the initial complex has been determined the 
iteraticn preceeds ty evaluating the objective functior J (x) 
at €ach vertex. The vertex at which the oktjective function 
has the largest value is then designated as the current 
worst vertex and this point is reflected througk the 


centroid of the remaining vertices establishing a new 
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W 
couplex. If this worst vertex is designated as x then the 
overreflecticn is accomplished according to the formula 
N CR wW 
xX = (1.0 + alpha) ex - alphaex 


where alpha 2 1.0 is the overreflection ccefficient, x 
Lepresents the centroid of the remaining vertices,and a is 
the new vertex. If the overreflection of the current werst 
vertex resuits ina violaticn cf the constraints or if the 
new point still has the worst value in the set of vertices 
then the feint is mcved halfway toward the centroid used in 
the overreflection process. This retraction is repeated (if 
necessary) until the overreflection coefficient is reduced 
to some stall value, delta, or a Suitable new vertex is 
established. If the overreflection coefficient is reduced 
to the value delta and the objective functicn value at this 
néw foint is still the worst in the set of vertices then the 
prcejected vertex is replaced by its original value and the 
second wcrst vertex will be overreflected instead. This 
prccess keeps the complex moving toward the wpinimun of the 
ckjective function provided the figure has ‘not ccllapsed 
mato its centroid. Once assured that the new ccmplex 
Satisfies all the constraints and results in an improvement 
of the okjective function, vertices of the complex are 
rejected and generated in a systematic manner aized at 
Binipization cf the objective function. Thus essentially 
the complex moves over the feasible region eventually 
straddlin¢c and collapsing upon the point that renders a 
Birimum value for the objective functicn. 


This frecess ccntinues and is eventually terminated when 


the complex shrinks to a predetermined acceptaktle small 
Value epsilcr given by: 
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k 
1 Cc v 2)0.5 
{ = et (eee ] < epsilcn 
1=1 


c 
Where J(x ) represents the value of the objective function 
V 
at the centrcid ard J(x ) is the value at the varicus 


vertices. 


As pcinted out by Box [4], the exact values cf the 
overreflecticn coefficient, alpha, and the exact numter of 
vertices, k, do not appear to be critical to the functioning 
of the algcrithm prcevided that they are greater than unity 
and nt1 respectively. The main concern here is that a value 
cof alpha larger than unity keeps the ccmplex from shrinking 
prematurely as it traverses the feasible region in search of 
the pinizum ckjective function value. Values of alpha = 1.3 
and k = zn as given in ref. 4 are used here since the exact 
values chosé€r are apparently not critical and good results 
have been reperted, by the originator of the algorithna, 
using these farticular values. For a detailed example of 
the iteraticn process and a more sophisticated explaration 
of the algcrithm the interested reader is directed to refs. 
2, 4, and 6. 
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A. DESCEIETICN OF ALGORITHS IKELEMENTATICN 


The frequency domain algorithm described in the latter 
part of Chapter II has been irplemented in program forn 
using FCRTEAN IV ccding and ccnfigured for use on the IBM 
360-67 ccrputer system available at the Naval Postgraduate 
Schcol Ccmputer Center. With the exception of the graphing 
sukroutine tke remainder of the program is self ccntained 
and may be used on any computer system having a standard 
FORTRAN IV compiler. Without the overhead routines 
reguired of the graphing subroutine, the basic prcgram 
reguires apfroximately 180k bytes of ccre memory, and 
including tke core requirements of the flotting routines 
roughly 210k ktytes of memory are reguired ky the prcgram. 
While admitedly less core would be necessary if such 
techniques as array overlaying were employed when 
prcgramming the algorithm, it was felt that for sinuplicity 
in following the flow of the program by future users these 
technigues were best left unused. 


The ccmputer aided linear time invariant compensator 
Ooptimizaticn pregram (CALICO) is essentially composed of 
fcur main farts with numerous additional small sukroutines 
functioning in supportive roles. Since often a prograe user 
May wish tc perform certain modifications tc programs of 
this nature in order that they may better serve his 
particular purpose, Appendix A provides a descripticn cf the 
SUffortive sukroutines used within the frogram to facilitate 
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acccmplishment of any program mcdificaticns that may be 
desired. A block diagram of the major components cf the 
overall program are shown in figure IV-1. The main program 
serves frimarily tc input data describing the system under 
consideraticr and to output the values of the roots and the 
gain of the compensator that will best achieve the desired 
response in accordance with the cost function selected. 
During program execution, the main frogram also serves to 
@Girect the seguence of execution of the other major sections 
of the CALICC algcrithno. 


Since the plant and compensator are assumed to be linear 
time-invariant, they may be completely described by transfer 
functions in the independent Laplace variable s. Once the 
flant and assumed f¢rm of the series compensation are input 
alcng witk the profiles forming the desired phase and 
magnitude respcnse program control Froceeds tc the 
Einimization subroutine, BOXPLX. From the point of view of 
computer aided ccntrol system design, asscciated witk this 
surroutine are two key function sSub-programs, FE(X) and 
KE (X). The function Sub-program FE(X) performs the 
calculaticn cf the system open loop frequency response, at 
the discrete frequencies provided in the input data, and 
calculates tke value of the cost function FEased upon the 
difference tetween the actual and desired frequency response 
at the discrete frequency values. Six separate cost 
functions, frem which the program user may select cne, are 
incerpecrated into this function sub-programp. The value of 
the cost function is used by the ninimization routine to 
vary the ccmpensatcr parameters in order to minimize this 
quantity. The function sub-program KE(X) checks fcr any 
Violaticn cf the implicit constraints, which in this 
particular case are concerned mainly with the presence of 
right half plane rcots. The iaplicit constraints used, also 
impose limits on the real and imaginary parts cf the 
huperator anc deno#ginator polynomials evaluated at any 


a0 


MAIN PROGRAM 


Output optimized results; call simulation, 
Stebilvty, and mlotting nemtines - 


SUBROUTINE BOXPLX 
Systematically varies the compensator 
parameters to minimize the cost function J. If 
the number of trials is exceedid the best 
solution achieved to that point is selected as 
He SOlutcion. ££ tne cost function remains 
unchanged for a soecific number of trials the 


minimum is considered as having been achieved. 


PUICTION FE PUNCT Lene KE 
Calculates the fre- Gnecks £Or samgnt 

yuency response of nNalf olane roots 
the open loop system ac Vlolacionsor 
at the discrete val- Cuber implTert 
ues svecified in the constraints. 
miouc. and forms the 
difference between 
gesiredsand actual 
response making up 
ice COStmeunction. 


SUBROUTINE SIMULA; Calculates the open loop | 
frecuency response using ontimized parameters. ‘< 


SUBROUTINE GRAPH; Plots the desired graphs 
Svcecified bv the user. 


Program Block Diagram 


Figure IV-l 
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particular frequency. These latter limits are a direct 
result of the finite computer word size and the particular 
methcd tsed to represent numerical values withir the 
machine's mencry. Bcth of these factors have an effect in 
deterpining the s2aximum and minimum size numerical value 
that may ke represented within any fparticular comfuter. 
While the ccefficient values determined by the minimization 
routine will nct exceed these waximnum or pinimum values, the 
real and isaginary parts of the numerator and denominator 
polynomials evaluated at a particular frequency are squared 
during the frocess cf determining the tsagnitude cf the 
frequency respense. Without the implicit ccnstraints these 
squared values can exceed the size of the numerical 
capability of the sachine, particularly in the case cf high 
crder systems, causing premature termination of the 
algcrithm ané hence failure in obtaining a solution. Cnce 
tke miniwizaticn subroutine determines the ccefficients that 
result in a ninimum value of the cost function cr the 
Ginimization process is terminated for other reascns, 
prcgram ccntrcl is transfered back to the main fpregran. 
Here the reots and gain of the compensator which resulted in 
a GBinimum cest function are calculated and cutput. The cpen 
locp freguercy response of the sytem is then calctlated 
using a ucre densely distributed number of frequency values 
than would normally be selected by the designer fcr use in 
the minimizaticn process. This is done largely to show the 
existence cf any resonant feaks which might otherwise be 
straddled and not apparent if only the discrete frequencies 
selected Fy the designer were used in the plotting routines. 
The stability cf the entire system is then checked by adding 
the open lccp numerator and denominatcr pelynomials and 
applying the Rcuth stability critericn tc the resulting 
FClynomial which is the system characteristic equation. 
Firally tke results at the discrete frequencies specified by 
the designer and the "continuous" curve generated Ly the 


SIKFULA surkrcttine are output cn Bode and gain versus _ fhase 
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picts. In addition, plots of the difference between the 
desired and actual results at the discrete frequencies 
specified by the designer are output in crder to aid the 
designer in sgaking any modifications to the form of the 
cc&rpensator that may be necessary in order tc better satisfy 


the system specifications. 


As stated in the previous paragraph six separate cost 
functions have been incorporated into the furction 
Suk-program FE(X). The first three of these are cf the 
€XLror squared type and the latter three are cf the aksclute 
error type. These six cost functicns represent various 
ccmbinaticns cf one minus the actual response divided Ey the 
desired respense at the discrete frequencies under 
consideraticr. These pre-defined cost functions are 
prcevided as a matter of convenience in using the pregram and 
certainly de nct prevent the user from defining his cwr cost 
function prcevided this is placed within the FE(X) function 
suk-program. The first cost function combines the deviation 
of koth tke desired magnitude and phase resfonses frem the 
actual magnitude and phase responses of the ofen Ilcop 
system. That is the first cost function is cf the fors 


> a ] (» “ee Y : (1.0 - me 
DMAG (a) DPHASE (wo) 
where AAG and APHASE are the actual magnitude and phase 
response cf the cpen loop system and DMAG and DPHASE are the 
desired sagritude and phase respectively. These quantities 
are evaluated in the cost function and summed cver the 
discrete frequency values in the specified range. The 
second ccst functicn which the user may chocse, uses solely 
the deviaticn of the actual magnitude response desired 
Magnitude response evaluated at the discrete frequency 
values. Thus this cost function is of the general fers 
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AMAG (w) : 
J (w) -) (1. see) 
DMAG (w) 


If a winisum phase system is being considered the magritude 
and phase flcts as a function of frequency are uniquely 
related ty the Bode theorems {3]. That is, determining the 
Magnitude ratio cf a system over a given frequency range 
will alec result in a unigue phase relaticnship for the 
systen. Thes if the systen design being considered is of 
the mininum phase type the use of the type two cast function 
described skculd result in both the magnitude and phase 
specificaticrs being satisfied. The third pre-defined cost 
function uses tke deviation of the actual open lcop fhase 
response fcrs the desired phase response in a functicn of 
the form | 


J (#) > (1.0 ae 
DPHASE (w) 
While for sinisum fhase systems the magnitude and phase 
responses are wuniguely related, it should te noted that the 
fixed lcss cr gain of the system cannot be determined from 
the phase characteristics alcne. Consequently the use of 
the type three cost function may result in a constant gain 
errer fcr tte magnitude response of the system. This error 
however is readily observable from the output produced by 
the program and correction of this may be accomplisted by 
adjustment of the compensator gain value accordingly. The 
cther three pre-defined cost functions are of the same 
general fcrgz as these just presented with the exception that 
the squared differences are replaced by the absolute value 
of the differences. The general rationale in selecting this 
latter type cf cost function was that on occasicn the 
Squared type cost function will have a tendency to flatten 
hear the gininmum. With the absolute value ccst functicn the 


user may avcid the excessive time required in findinc the 


ay] 


Binimum valtes wten the case arises where the squared cost 


furcticn Leccnes flat. 


In utilizing this design program the user should alsc be 
aware that the need fer frequency scaling cf the problem may 
arise. This is particularly true when high order systems 
involving high frequencies are being considered. The need 
for freguency scaling is primarily a result of the finite 
word size available in the computer system. It must Le kept 
in mind thet the magnitude response of the system is 
computed as the ratio of two ccmplex numbers and while the 
resulting macrnitude may be relatively small the value cf the 
Bagnitudes of the individual complex numbers forming this 
ratio may te large. The situaticn is further complicated by 
the fact that in ccmputing the magnitudes of the ccamplex 
numkers invclved the real and imaginary parts must be 
eguared, thus effectively reducing the maximum power of the 
numbers tshich can be allowed by a factor of two. If the 
magnitudes of the resulting numkers kecome too Jlarce an 
errcr messace will ke printed recommending that the problen 
be freguercy scaled. 


B. DATA INECT ANI OUTPUT 


In crder to sigplify use cf the program, an attemft has 
been made tc maintain the input information necessary for 
executicn cf the algorithm to a minimum, and at the same 
time preserve a certain amount of flexibility withir the 
Frogram which may be controlled by the user. In so doing 
scme of the program control variables for various cftions 
have been set irternal to the program and may not be 
Manipulated ty the user in the input data deck. These 
internally set variables do not detract from the general 


functionirg cf the program, but mention of this is made to 
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bring it tc the user's attention in the event that 
specialized frceblenms might more effectively be investigated 
by resettirg scme of these normally preset prcgranm 
variables. Many of the preset items are concerned with 
centrcl cf the tinimizaticn routine. Such items as the 
upper and lcwer bounds of the search area, the integer 
prcgramming cption, and the seed for the randcm number 
generator are set to specific values within the main frogran 
just pricr tc calling the minimization subroutine. For 
example fcr the aginimum phase case, the lcwer and upper 
bcunds cn the search areas of the compensator numeratcr and 
dercminatcr fpelynomial coefficients are set at 0.C and 


6 
1.€x10 respectively. For the non-minigum phase case the 
beunds on the seé€arch areas cf the hnumeratcr pclynonial 
6 6 
coefficients cnly are extended to -1.6x10 and 1.6x1C and 


the implicit ccnstraints are adjusted to allcw fright half 
Flane rccts for the numerator polynomial. The integer 
prcgramming cfpticn for BOXPLX is also preset:'to a value of 
z€xro indicating that values for the pcliynomial coefficients 
need not ke integers. In general, selection of the irteger 
prcgragming cpticn results in considerably more execution 
time fcr the program and for the purpose cf varyinc the 
compensatcr farameters the selection of this option precvided 
ncthing tc erhance the soluticn of problems. 


The input data reguired for program executicr and 
control may te divided into general areas as follows: 1.) a 
title card; 2.) a series of cards giving the gain and 
ccefficierts cf the plant transfer function; 3.) a series 
of cards giving the assumed gain and coefficient values of 
the compensatcr transfer function; 4.) a control card 
describinc the flots desired, the type of cost functicn to 
Fe used, whether the system contains a zero crder hold, and 


tke number and range of the discrete frequency points being 
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considered; £.) cards containing the discrete frequency 
values, the desired magnitude, and the desired phase; 6.) 
finally a card describing the number of iterative trials to 
be allowed in the minimization and whether any cutfut is 
desired frem the minimization routine in order to keep track 
of the sganner in which the coefficients are being varied. 
The specific input formats of the above necessary data are 
described ir detail in the following faragraphs. These 
paragraphs are numkered so as to indicate into which cf the 
six general areas a particular data card belongs. This 
descripticn ccupled with the program listing given in 
Appendix 2 shculd provide sufficient informaticn for setting 
up the data deck necessary for execution cf the frecgran. 
Figure IV¥7r2 illustrates the data deck used in executing the 
first example prcblem presented in section C of this 
chapter. In this figure each line of data represerts an 
incividual data card. The column numbers for the data cards 
are alsc shcwn abcve and below the data deck secticn cf the 
figure. 


RS — << 


dq. Title card. 


This must be the first card of the data deck and 
Frovides a s#eans of problem identification if several 
s€parate pretlems are crcun under one jok name. Columns 1 
through 4E may be used to input whatever alphax~<numeric 
characters the user desires for identifying any particular 
Froklem. This title will also he printed cn the grafhical 
output reguested by the user. 


The seccrd card contains the gain of the plant in 
cclumns 1 through 10. If the gain is not explicitly present 
in the plant transfer function then the value input cn this 
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card shculd ke unity. The format used tc read the plant 
gain is ar £10.00 type format. Thus, the gain may be input 
in either floating point or scientific notaticn type 


formats. 


This card specifies whether the numeratcr of the plant 
transfer functicn is to be read in factered cr pclyncmial 
form and what the crder of the numerator is. A P or an F in 
cclumn 1 indicates whether the plant numerator data is tc be 
read in pclyrcmial or factored form respectively. Columns 2 
and 3 cortain the crder of the numerator. Thus an input of 
the form FC3 starting in column 1 indicates that the 
EClyncmial is of third order and is to be read in factored 
fora. 


The fclicwing card(s) contain either the pclyromial 
coefficients cr the factors of the plant transfer function 
nuwerator, depending on which form has teen chcsen for 
put . If the pclynomial form of input has been selected 
then the gEclyncmial coefficients are input in ascending 
Powers of ¢ using an 8£E10.0 format. If more than eight 
coefficierts are needed they are simply continued on 
successive cards. The coefficient of the highest order tern 
is assumed tc be nermalized to unity and is not read in 
explicity. Thus the number of coefficients actually read 
ccrresponds tc the order of the numerator folynomial. rf 
the factcred form of input is chosen then each factor is 
read in cn a separate card with columns 1 through 10 being 
used for the real part of the factor and cclumns 11 thrcugh 
20 for the imaginary part. Again the format used is £10.0, 
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thus allcwirtg the factors to be input in either flcating 
peint or scientific notation forn. Since complex factors 
will occur in ccnjugate pairs the values cf these factors 
need only re input cnce. Thus the number cf cards used will 
vary deperding ufon how many purely real and how many 


corplex factcrs are involved. 


_——_ = == ee a a _— Se ee 


Again tke fcrm of input is specified as in card number 
22 along with the order of the plant transfer furction 
dercminatcr. 


lant transfer function derominatcr. 


— _—= ae ae oe =a a ae So = Se ee ee ee 


The fcrmat used to input the denominator is identical to 
that used fcr input of the numerator, contingent upfen the 
desired fcrm specified for the input. That is to say, it is 
possitle tec input the numerator in polyncnial form and the 
dencminatcr in factored form, or vice versa, since each tine 
a folyncsial is input as data the input form is uniquely 
specified. 


Folicwing the input of the information about the plant 
the descrirfticn cf the form cf the compensator must be 
Supplied. This is begun by inputing an assumed comrersator 
gain using tke same format as was used to specify the plant 
gain. Cue to the manner in which the compensator transfer 
function nugeéerator coefficients are handled withir the 
Frcgram, if the assumed compensator gain is set to zerc, for 
lack cf any fetter value, whis will result in the starting 
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guess of all the numerator coefficients being zero. fhis 
will not cause the program to terminate prematurely, but the 
initial ouess of zero for all the coefficients may not be 
what is desired by tke user. If doubt exists as to what 
value gain shculd te assumed, a value of 1 is suggested as a 


convinient arkitrary starting point. 


The methcd for specifying the form of the input fcr the 
ccmpensatcr initial numerator parameters is identical as 
that used in specifying the plant numerator. A P indicates 


peclynomial fcrm and an F indicates factored forn. 


In crder to maintain continuity in the input fcrmat the 
scheme and fcrmat used to input the compensator numerator 
parameters is identical to that used in specifying the flant 
transfer functicr numerator. 


As in frevious cases either polyncemial cr factcrec¢ form 
is specified along with the crder of the compersator 
dercminatcr. 


Once again the same format as used to input the flant 
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éata is used in order to maintain the ccntinuity cf the 


Moeut forsatirg. 


fhe precram control card is used tc input various 
quantities ccntrclling both the range and number of discrete 
peints used in the minimization of the ccst functior, the 
output desired frcw the program, and also whether a zero 
order hcld is present in the error channel of the system. 
As many as twelve separate quantities may be input on this 
sirgle ccntrel card in order to direct executicn cf the 
program. Cclumns 1 through 10 specify the minimum frequency 
(WMIN), in radians, of the total range Lteing consicered. 
Cclumns 11 through 20 specify the maximum value (WMAX), in 
radians, cf the range of freguencies being considered. FPoth 
of these values are read using a 2E10.0 format. Fclicwing 
the waximun and minimum frequency values, up to ten input 
quantities are read using a 1013 format. Columns 21 through 
23 specify the numker of discrete frequency foints (NOMEG) 
over the given range at which the open loop frequency 
response is to be determined and the cost furction 
evaluated. The next input quantity (KNOW), in columns 24 
through 26 specifies whether the discrete frequencies of 
interest are to be read from data cards or automatically 
computed. Ncrmally this quantity is set to unity indicating 
that the discrete frequency points of interest will be read 
frem data cards. A value of zero indicates that the 
discrete frequency values are to be incremented linearly by 
the program starting at the minimum frequency value. A 
value of two indicates that the discrete frequency feints 
will be incremented logarithmically begining witk the 
tinimum fregvency value. Columns 27 thrceugh 29 indicate 
whether a Ecde flot of the results is desired. If this 
variakle (NECCDE) is set equal to one no Bode flot will he 
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output. A value of zero for NBODE will result ir roth 
Bagnitude ané phase flots. Columns 30 threugh 32 ccntrol 
the output cf the Nyquist plot. If NYQST is set to unity 
then the Nyguist plot will ke suppressed. A value cf zero 
for NYQST will result in a Nyquist plot output, however this 
is plotted crly over the range of frequencies between WHIN 
and WMAX. Cclumns 33 through 35 indicate the presence of a 
zero crder hcld (1Z0H) in the error channel of the unity 
feedback system. With IZOH set to zero the system is 
treated as keing of a continuous nature. With IZOE set 
egual tc cne a zero order hold is considered to be present 
in the error channel and the open loop magnitude and fphase 
computaticns are modified acccrdingly. Cclumns 36 through 
3& (NICHCL) specify whether a magnitude versus phase plot of 
the results are desired. As with the other flot opticns, a 
value of zere will result in a plot being preduced while a 
value of cre will suppress the magnitude versus phase flot 
as part cf tke output. Columns 39 through 41 contain the 
variable indicating whether non-minimum phase soluticns are 
tc be considered. With this variable (NMINFS) set equal to 
zerc only sinimum phase solutions wiil be allowed. A value 
of cne indicates that a non~minimum phase solution will be 
allowed provided that it is the one which rinimizes the cost 
functicn. Cclumns 42 through 44 indicate if the desired 
magnitude prcefile values tc be input are in decibels. If 
this variable (IDB) is zero the magnitude values to te read 
are disensicnless gain values, while a value cf one 
indicates that the magnitude values to be input will fre in 
decibels. cCclumns 4S through 47 contain the variable IPLCT. 
This guantity is used to choose between printer flots or 
Calccmp flcts for the graphical output. A value cf zero 
will result in printer plots and a value of unity will 
previde cCalccmp cutput of the graphs desired. It is 
reccmmended that the printer plot opticn be used in 
initially sclving a particular problem since this type of 
output recuires less turn around time than the Calcomp 
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output. Finally columns 48 through 50 specify the tyre of 
cest function (ICOST) the user wishes to chcose amoncg the 
six preprogrammed functions available in subroutine FE(X). 
The value specified should be between 1 and 6 correspcnding 
te the farticular cost function desired. If this value is 


zere the procram will default to the type cne cost function. 


Sampling period. 


This card must be inserted only if a zero order hcld is 
specified in the error signal channel. The sampling feriod 
T, in seccnds, is specified in cclumns 1 through 10. If the 
system btéing ccnsidered is of the continuous type then this 


card must ke critted. 


Next fcllcws a set of data cards specifying the values 
in radians fer second of the frequencies at which the 
desired gain and phase values of the open loop frequency 
response are specified and at which the cost functicn is to 
be calculated. These values are input in an 8£10.0 format 
using as many cards as necessary to speciry all the 
frequency values at which the program is to execute the 
alccrithg. 


This set cf data cards contains the desired ofen loop 
gain respconse at the discrete frequency values previously 
specified. MThe samze 8E10.0 format is used here as was the 
case in reading the discrete frequency values. The same 


nugker cf cards shculd be used as was used in specifyirg the 
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discrete frequency values. 


Here fcllcws a series of cards specifying the desired 
open locg fhase response, in degrees, at the discrete 


freguency values being considered. 


Here ktoth the number of trials to te allowed in the 
Zinimizaticn frocess and the print interval for diagrestic 
purpeses are specified. These quantities are input using a 
215 format. If the number of trials is zero then a default 
value of 2C0C is assumed. If the print interval is 
specified as zero then no output from the minimization 
routine will result. Caution should te exercised in 
Specifyinc the print interval in that if this interval is 
made too snall an excessive amount of printed output may 
result. The program user is referred to the comment section 
of subroutine BOXPLX in Appendix A for guidelines in 
selection of the print interval. 
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Cdeeeeday 


Figure IV-2 


C. SERIES CCMPENSATOR DESIGN EXAMPLE FROELEMS 


In this section a number of example problems are 
presented and discussed in order to illustrate the use cf 
the computer aided compensator design algcrithm. In the 
hcpe cf presenting a clear, concise illustration of the 
algoriths and its use the example problems begin with 
relatively simple straight forward text bock examples and 
gradually regress tc problems of a slightly more difficult 
nature. Alsc included are problems that are intended to 
illustrate scme limitations of the program in that these 
prctlems were not completely solved by the algorithm cr they 
may have reguired a certain amcurt of manipulation cf the 
input data in crder to obtain a soluticn. These latter 
types of freklems have been included since it is impcertant 
that the ¢esigner using any computer aided design algcrithm 
understand the algorithm's limitations and peculiarities as 
well as its ixtended capabilities. It is nct intended to 
gislead the reader by implying that the example frecblems 
presented here provide a comprehensive presentation of all 
Fossible Jlimsitaticns of this automated design prcgran. 
Rather ttey represent those limitations which were 
enccuntered ir the process of verifying the program cver a 
necessarily linited period of time. With this in mind the 
first example precblem which deals with a relatively simple 


lead compensatcr is presented. 


Consider the uncompensated system shown in figure 
av -J. 
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4 
s(st+l) (s+2) 


Design Example 1, Block Diagram 


Figure IV-3 


The asymptctic Bede diagram for the open loop system is 
shown in ir figure IV-3A. As can te seen frem the cpen loop 
Bode plot, even thcugh the system is stable it is lightly 
camped. A single section lead compensator is to be used to 
increase tke phase margin. A complete discussion of the 
reasons for selection of the lead compensatcr and a detailed 
discussicn of this particular example is presented Ly Thaler 
and Brown [18]. The transfer functicn fcrm of the single 
section lead ccmpensator consists of a first crder numerator 
and a first order denominator. A desired magnitude and 
phase prcfile, consisting of ten discrete frequency foints 
over the range from Q.2 to 20.0 radians, was selected to 
give a phase margin cf approximately 45 degrees in crder to 
increase tke damping cf the system. As an initial guess the 
ccipensator was assumed to have a gain of 1.0, a zero at the 
Origin ard a ¢cole at ~1.0. This informaticn was suprlied to 
the programs in the prescribed format and the oftimized 
results in the form of magnitude and phase plots are shown 
in figure IV-3£. Here, as thrcughout the remainder of the 
example prcklems, the diasond shaped symbcls represent the 
desired values specified by the program user, the X symkols 
represent the values computed by the program at the discrete 


50 


frequency values specified using the compensatcr farameters 
returned as the sclution from the minimizaticn routine, and 
the continuots curve represents the response computed ky the 
pregram again using the parameters returned as a solution 
but with a mcre dense set of frequency values used in crder 
to show any resonant peaks or anomalies in the magnitude and 
phase curves which might be difficult to visualize frcem a 
sparsely selected set of discrete values. As can be seen in 
figure IV-32 the results match the desired response well 
over the trarge cf interest. Figure IV-3C shows the same 
results using the magnitude versus phase plct to portray the 
open locp frequency response. The numerical values cf the 
ccmpensatcr fparameters returned from the program as a 
scluticn alcng with the starting values used are shcwn in 
meaqures J¥-3C, IV-sE, and IV-3F. A comparison with the 
scluticn values given in ref. 18 will show that they agree. 
Figure IV¥-3¢ shows plots of the differences between the 
actual ard desired open loop frequency response cf the 
resulting system at the discrete frequencies considered in 


tke minirizaticn process. 


This problem as descriked reguired approximately 500 
iterations sithirn the minimizaticn routine to achieve the 
sclution <shcewn. This however includes twe restarts to 
attempt insuring that a global minimum had keen achieved and 
that terminaticn was not due tc a local mininun. Frogz the 
begining tc the first restart required approximately 250 
trials with the cost function being reduced from a value of 
6.23 to a value of 0.00711 (recall that the minimun 
achievable urder ideal condition is zerc which due to 
machine errers alone in representing numbers would not 
likely ever ke achieved). The entire soluticn requireé 32.5 
seccnds cf CFU tine, a good 10 seconds of which was used in 
the prodtcticr of Calcomp plots. These numbers ccncerning 
the time required for soluticn are presented for what they 
are wortk in that extrapolaticn of these numbers tc more 
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ccuplicated problems is not a straight fcrward task. As 
menticned in Chapt III the minimization routine does [frecome 
inefficient when the number of variable parameters is 
allewed tc kteccme large. It was felt that a detailed 
investigaticn of the efficiency and time required for 
sclution cf froblezs as a function of the number of variable 
parameters and the number of discrete frequency values tc be 
censidered was beycnd the scope of the kasic developmert of 
the progragr. However, the amount of time required for 
sclution cf scme of the example problems that follow will be 
given in crder to give the user an "intuitive feel" fcr the 
approximate amcunt of time that may be need to kandle 


certain tyres cf problems. 


In crder to illustrate that the starting foint 
chcsen shculd have reasonable values, but that the exact 
values chosen should not effect the final results 
significartly, the same problem was run a second time with 
the pole and zero of the compensatcr both chcsen tc fe at 
the crigin. Again an assumed initial gain cf one was chcsen 
for the ccmpensator. The graphical results cf the solution 
are shown in figures IV-3H through IV-3J. The numerical 
results returned using these new starting values are not 
shewn here since they were identical to within three decimal 


places with those values shown previously. 
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2. Ccmrensator Lesign, Examp 


lo 
Iho 


q{kis example problem (discussed in detail in ref.10) 
has been chcsen nct only to familiarize the reader with the 
use cf the program, but also to illustrate one of the 
shortccmings cf the design prcgram. In particular, while 
the program was successful in finding the correct sclution, 
difficulties were encountered in the root finding 
subroutines. In this particular problem, during the final 
phases cf execution the roots cf a fourth crder pclynonmial 
must ke fcuré in order to calculate the fhase response using 
sukroutine SIZFULA tvhich does these calculaticns using a more 
densely spaced set cf discrete frequencies than ir the 
pinimizaticn routine. Since numerical root finding routines 
have difficulty in finding roots for fourth crder 
Eclynomials, the lack of good convergence for the rocts of 
the system in Ecth rcot finding subroutines caused the phase 
calculaticns normally done by subroutine SIMULA tc be 
abcrted. It should be emphasized that failure of the root 
finding subrcttines in this case only effects the generation 
of the "ccntinuous" phase response and in nce way created any 
prcklems fer the sminimizaticn part of the prcgram as 
evidenced Ly the sresulting correct values for the 
ccmpensatcr farameters. The original unccmopensated systen 
being ccrsidered is shown in figure IV-4. A Bode diagram of 
the unccupensated system (figure IV-4A) reveals that the 
system as it exists is unstable. If stability were the only 
consideraticr, a simple attenuation of the plant gain would 
be sufficient tc provide adeguate phase margin tc insure 
this. fewever, as with many series compensation frotlems, 
here we must effectively strike a compromise between 
effectiveness cf ccntrol (generally asscciated with a higher 
gain) and stakility of the system (generally associated with 


a lower gain). As discussed in ref. 10a simple gain 
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attenuaticn wculd result in stability being achieved, but 
also has ar undesirable effect upen the velocity error 
constant. Tkus as discussed in detail in ref. 10, a lag 
compensation retwcrk appears as the best candidate to 


acccmplish tke job. 


5°+20s°+100s 


Design Example 2, Block Diagram 


Figure IV-4 


In this particular problem twenty discrete frequency 
values over the range from 0.005 to 50.0 radians were 
selected ard the desired magnitude and phase values 
specified at these frequencies. As in the frevious eéxample 
proklem the form cf the compensator is selected as having a 
first crder numerator and first order denominator, that is 
in this case a sincle section lag compensatcr is to ke used. 
A type one fcrm of the cost function was selected for the 
Binimizaticn process. The magnitude and phase profiles of 
the desired response along with the resulting response with 
the compensator included in the cpen loop system are shown 
in figures IV-4B, IV-4C, and IV-4D. As mentioned, the 
continuous curve for the phase response was not computed 
because cf the failure of the reot finding recutines to find 
the cpen lccp system roots. It can be seen in figure IV-4C 


that the phase response does agree well at the discrete 
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frequencies specified. Since we are dealing with a pininuo 
phase system and figure IV-4B shows that there exists no 
resonant freaks in the magnitude curve the user may conclude 
that the phase response will also Fe well tkehaved fcr _ the 
various frequency walues ketween the explicitly specified 
discrete values felotted. This case should SEIVE to 
illustrate that the program reguires that the user still 
prcvide ar interpretation of the results. The coaputer 
output cf figures IV-4E, IV-4F, and IV-4G shows that the 
ccgapensatcr farameters are in agreement with the freblen 
scluticn as resented in ref. 10 if the reader wishes to 
make the ccmparison. Figures Iv-4H and IV-4I show flots of 
the differences between the actual and desired magnitude and 
phase at tke discrete freguencies selected for the 
Binimizaticn frocess. 
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Anctker lag type compensation protlem is fresented 
in this example. This time, however, the existence of a 
fourth crder peclynomial for e¢ither the cpen loop systen 
transfer functior rumerator cr denominator has furfgosely 
keen avcided tc insure that the program will run te 
cctpleticn. The uncompensated system blcck diagram is shown 
in figure IV-5. 
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Design Example 3, Block Diagram 


Figure IV-5 


A Bcde diagras for the uncompensated system along with a 
detailed discussicn of the problem may be fcund in ref. 17. 
The unccapfensated system does not possess an adequate fhase 
Bargin tc irsure froper system performance. Thus, a single 
secticn lag ccmpensator is to be used in ecrder to achieve 
the desired cpen loop frequency response. FCesired gain and 
Fhase prcefiles are selected for the ccmpensated ofen loorp 
frequency response using 20 discrete frequency points over 
the range frem 0.01 tomes00. 0) Ladians. The prcfiles 
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selected alcng with the results achieved are shcwn in 
magures I¥-SA, IV-5B, and IV-5C. The nuperical values 
assumed as initial estimates of the compensator parameters 
and the numerical results of the output are shown in figures 
IV-SD, IV-SE, and IV-5F. Plcts of the differences Letween 
the actual magnitude and phase, and the desired magnitude 
and phase are shown in figures IV-5G and IV-5H. In this 
particular prcklem approximately 1200 iterations of the 
tirimizaticn routine and 1 minute and 25 seccnds were 
required for sclution of the _ problen. In order to 
illustrate the idea that for minimum fhase systems the 
ccmpensatcr parameters may be determined to within orly a 
fixed lcss cr gain when using only the phase profile, this 
same problem was solved again using the type 3 cost furction 
which ccnsiders only the difference between the desired and 
actual phase at tke discrete frequencies. The results are 
shcwn in figures IV-5I, IV-5J, and IV-5K. As can be seen 
frcm the magritude plot in figure IV-5I, there is a constant 
gain errcr in the resulting magnitude curve. The numerical 
values assumed at the start and those returned after the 
tinimizaticn cf the cost function are shown in the couputer 
output of figures IV-5L, IV-5M, and IV-5N. The resultant 
Magnitude response is higher than desired. The necessary 
cerrecticn tkat must be applied to obtain the desired 
response may ke read directly from figure IV-50, which shows 
the difference between the actual and desired magnitude 
response. That is, as can be seen from figure IV-5C, the 
Cesultant tagnitude of the ofen loop systen after 
ccupensaticn is approximately 36.8 db higher than what is 

desired. Since the program does not alter any of the plant 
Farameters this peans that the compensator gain value 
returned frcem the program when the type 3 ccst functioE was 
used must be decreased by approximately a factor of 69.2. A 
cc@fariscn of the compensator gain values shcewn in figures 
IV-S5SF and IV-5N will show that they due indeed differ by 
this value. Also the astute cbserver May nctice that there 
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is a slight difference in the compensatcr pole value 
returned as the sclution when the problem was solved using 
the type 3 cost furction. While the difference in valves is 
small it is suspected that in this particular case it is due 
to the fact that only a finite frequency range has been 
considered ard the lower frequencies have been "slighted" 
scuewhat in that the nonlinear phase curve has nct yet 
flattened cut at 0.1 radians. This is also the same 
situation tkat existed the first time the prcgram was 
executed, but in this case the magnitude prcfile was also 
taken into account in the cost function and any deviations 
cf the resultant magnitude values dominated the cost 
functicn value. In other words, it appears that ir this 
particular protlem the cost function is more sensitive to 
variations in the magnitude values than in the variaticns in 


the phase values. 
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Figure IV-5N 
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In this example a slightly more complicated problem is 
presented, in that the final soluticn requires a double 
section ccmpensator of the fcrm of a notch filter. A _ block 
diagram cf the unccmgensated system is shown in figure IV-6. 
As discussed Ey Thaler and Frown [18], from whom this 
example was taken, the uncompensated system is unstable and 
the Bcde diacram of the uncompensated system shown in figure 
IV-6A indicates a phase margin of approximately -30 decrees. 
With such a large negative phase Margin £Cr the 
uncompensated system and the rapid drop in slcpe cf roth the 
gain and fghase curves the designer might suspect that 
compensation using only a single secticn ccmfensator will be 
difficult to accomplish. Again desired magnitude and fhase 
prefiles are selected and initially a single section of 
ccmpensaticn is assumed to observe how close to the desired 
response tke system will perform. As can ke seen in the 
resultant craphical putput of figures IV-6A, IV-6B, and 
Iv-€C, the single section compensator fails to meet the 
reguired frecuency specifications. In fact, while it is 
Ecssible to stabilize ths system with a single section 
compensatcr the bandwidth will be excessively large and the 
desired Bagnitude profile prevents the program fren 
acccmplishing this. The numerical results returned frem the 
Single secticr ccmpensator run are shown in figures IV-6D, 
IV-6E, and IV-6F. As can be seen in figure IV-6F, the Routh 
test of tke characteristic equation also indicates systen 
instability which is to be expected. The magnitude and 
Phase difference curves of figures IV-6G and IV-6H, which 
indicate tte difference between the specified magnitude and 
phase prcfiles and the values actually achieved, indicate 
that additicnal compensation may be needed at the higher 
fréeguency ranges. Thus a double secticn compensatcr was 
assumed and the resudts achieved are illustrated in figures 


94 


IV-61, IV¥-6J, and IV-6K. The program is able to satisfy the 
specificaticns with a double section of compensation and the 
parameters required to accomplish this are shown in figures 
IV-6L, IV-6¥%, and IV-6N. The differences hbetweer the 
desired and specified magnitude and phase values fcr the 
double secticn cf compensation inserted in the system are 
shcewn in figtres IV-60 and IV-6P. 
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Design Example 4, Block Diagram 


Figure IV-6 
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this final compensator design example represents a 
prcblem cf a slightly different nature than those fresented 
to this fecint. The servo system under consideration here 
bas a zero order held present in the error channel ard in 
addition, due to the large coefficient values that result in 
the plant transfer function, the problem requires frequency 
scaling in crder to insure that the numerical limitaticns of 
the particular computer that was being used would nct be 
exceeded. {he original single loop system is shown in the 
block diagrams cf figure IV-7. Series compensation of this 
system was desired in ocrder to meet the fcllowing 
specificaticns; 1.) the open loop magnitude at S07 
radians/seccrd 2 30db, 2.) the gain crosscver frequency 2 
18E5 radiansysecond, and 3.) the phase targin at gain 
crcssover 2 45 degrees. These specifications were to re met 
using only a series compensator and the parameters of that 
part of tke system shown in figure IV-7, other than the 
compensatcr, waS tc remain unchanged. An initial attempt at 
using the CALICO program to design a ccmpensator resulted in 
a warning of excessive number size during computaticr and 
reccmmended frequency scaling of the problem. Thus a scale 
factor of 10C0 was selected as a convenient value and the 
entire system, including the sampling frequency of the zero 
order hold, sas scaled down in frequency by the valve of 
this scale factor. Using Laplace transfcrm notaticn, this 
Was acccaplished by a straight-forward substituticn of 
S = 1000S, where S represents the new scaled freguency 
values. A bleck diagram shcwing the values of the systen 
Parameters after scaling is given in figure IV-7A. When 
this type of scaling is performed the resultant numerical 
values cf the program need only be multiplied fry the 
appropriate ccenstant value in order to be corrected for use 
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in the unscaled systen. 


Criginally, work on this particular system indicated 
that a fifth order compensator could provide the necessary 
phase margin and gain crossover frequency, rut that the 30 
db specification at 377 radians/second was not satisfied. 
Thus, with this information as a starting point the 
Magnitude frofile cf the system with the fifth order 
ccmpensatcr in the loop was adjusted at the .low frequency 
end to satisfy the 30 db requirement at the scaled frequency 
value of C.377 radiansysecond. The phase prcfile was left 
unchanged. The results of the program can be seen in 
figures IV-7E through IV-7I. As can be seen frog the 
graphical cutput the desired specifications have teen 
satisfied ky the resulting compensated systeg. 
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Figure IV-7D 
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A. GENEBAL CLISCUSSION 


In tke frocess of running various test problems to 
exercise and verify the algorithm presented in the preceding 
chapters, it was found that in additicn to its use as a 
ccmfensatcr design program, the basic method could also be 
applied tc the synthesis of transfer functions from 
frequency respcnse data. That is, if the magnitude and 
phase measurements of some system are ocktained it is 
possible tc vse the program (CALICO) in order to oftain an 
approximate linear analytical expression for the transfer 
furcticn cf the system from which the measurements were 
oktained. This type of problem is guite cften encountered 
in. the mcdeling precess of dynamic systems. That is, the 
engineer may have measured input-output data of the system 
in the fcerm cf a frequency response and is confronted with 
oktaining a transfer function which will mcdel this system. 
This may ke the result of the process being toc complex to 
be analyzed and the pertinent equations written on the basis 
of physical laws or the governing physical laws may nct be 
completely understood in the case of systems involving new 
technologies or reasearch effcrts. Quite often the 
prccedure tcllowed in synthesizing or determining an 
approximate Jinear transfer functicn from frequency response 
data is te ccnstruct a Bode plot for the measured data and 
to fit the measured characteristic with straight line 


asymptotic approximations where possible and quadratic 
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factors where sharp peaks in the measured data occur. As in 
ccppensatcr design by classical frequency domain techniques, 
this precedure for synthesizing transfer functions can 
invclve a ccnsiderable amount of time consuming and tedious 
"cut and try" cn the part of the engineer. 


The algcrithm previously presented is so structured that 
the use cf tke program for this type of transfer furction 
synthesis is possible with no modifications to the coding 
itself. Since the algorithm was not originally developed 
for the specific purpose of transfer function synthesis 
there are several precautions that must be taken in crder to 
insure that the transfer function obtained does indeed 
adequately represent the system from which the measurements 
were oktained. In particular when using the program in 
crder tc synthesize a transfer functicn, a sefarate 
Simulaticn rcegram for time domain response should be used 
to verify the transfer function obtained. MThis procedure of 
ctecking the time domain response should be followed fcr any 
transfer function cbtained from frequency response data. 
The primary ccncern here is that it is possible to 
accurately represent the frequency response cf some systems, 
over a finite freguency range, by a transfer functicn of 
higher order than the true transfer functicn, but in some 
cases tke dcnrinant response to standard tise domain inputs 
will be ccnsiderably different than that of the actual 
systen. Alsc, while there is a direct relationshif between 
the time dcme2in and freguency domain responses througk the 
inverse Fourier transforgo integral, this integration, 
strictly speaking, is taken over the entire frequency range 
frcm -~0C tc +0, In practical applications it is ckviously 
necessary to limit the range of frequencies to be measured 
to some firite region. If, however, this regicn is not 
sufficiently large for a particular system the resulting 
transfer function obtained frem the frequency respcnse data 
Ray not acctrately describe the systen. fiwethis is the 
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case, it should become readily apparent in a time dcmain 
simulaticn cf the equation obtained from the frequency 
response data. {The possible pitfalls in the use of this 
algcrithm tc accomplish transfer function synthesis, such as 
the one just mentioned, will be discussed later when some 


example protlems will also be presented. 


A linear time invariant system may Le expressed as a 
ratio cf two frequency dependent pclyncmials of the forn 


nD n= 
as +a s + eeee Sams * sa N(s) 
n n=1 1 0 
BH (s) 3 ra = = 
bes te 5 s + ceo t > S “eb CD(s) 
m n> 1 1 0 


where s is the frequency dependent ccmplex variable cf the 
Laplace transfcrm. Thus the value of this transfer function 
is dependert upon the value of s and the unknecwn 
coefficients cf the numerator and denominatcr. Note that 
this is analcgous to the situation that occurs in the Cceésign 
of a cascade ccmpensator for a unity feedback system with 
the plant represented by a unity gain block. fhus fcr the 
purposes cf using the program, if we imagine the measured 
frequency response data to represent some desired system's 
open loop frequency response profile, and represent the 
Flant of that system by a constant gain of unity, then the 
transfer function H(s) to ke synthesized may ke considered 
the same as a cascade compensatcr in the system. Thus by 
representing the problem in this manner the CALICO program 
may be used to approximate the transfer function cf the 
actual systesz Ey varying the parameters of this pseudo 
ccmpensator (actually the transfer functicn of the systen) 
in crder tc minimize the cost function. The cost function, 
bere represents the error LELetween the mneasured frequency 
response data (entered as the desired rofile) and the 


fréeguency resfense simulation calculated by the prcgran. 
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Thus in effect the program minimizes the errcr between the 
transfer function approximating the system and the measured 
freguency respense of the system at the discrete frequency 
points which form the profile infut into the progran. 


As is the case when using the program tc design 
ccupensatcrs the user must assume a specific order fcr the 
numerator and dencminator representing, in this case, the 
transfer function to be synthesized. While unfortunately 
the programs will not vary the order of the numerator and 
denowinatcr gclynomials in order to obtain the best match to 
the measured frequency response, the user can employ the 
measured Gata in order to make an initial guess, as it were, 
of the form cf the transfer function. After this, one may 
scrutinize the results from the program, which grafts the 
measured and simulated responses, to make a decision if the 
fcrm of tke transfer function should be adjusted in crder to 
more closely represent the system. Here alsc, a time dcmain 
Sigulaticn sight ke extrerely useful in guiding decisicns 
regarding the particular form cr modification that should he 
tried in crder to obtain a "good" model for the system under 
investigaticn. The fact that the user selects the specific 
crders cf the transfer function numerator and dencminator 
and the minisizaticn routine varies the parameters to cbtain 
a Binimus ccst function while leaving the crders unchanged 
presents the potential of obtaining a low crder model of a 
higher crder system. While for linear systems, with known 
transfer functions, there are numerous analytical techniques 
available for accoaplishing this, the capability exists here 
of essentially oktaining a reduced order transfer function 
without a friori knowledge of the higher crder transfer 
function describing the system. The only data required by 
the program is the frequency response of the higher order 
System and ncthing concerning the parameters of the high 
crder transfer function is required. 
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Since scme transfer functions may contain roots of the 
nuserater feclyncmial in the right half of the s~plane 
(ncnminisum phase systems) this type of transfer function 
should rect ke overlooked when attempting tc model a systen. 
fo this erd a nonminimum phase cption may be chosen ir the 
CALICO cceputer algorithm. When this opticn is selected the 
search area fer the numerator coefficients is expanded to 
include Eoth positjve and negative values and the implicit 
censtraints requiring left half plane rcots for the 
nuperator efeclyncmial are eliminated. Selection of this 
cption dces nct guarantee that a nonminimum fphase transfer 
function will result from the CALICO program but it does 
allow this type of transfer function to be considered as a 
feasible sclution cf the problen. 


B. TRANSFER FUNCTION SYNTHESIS: EXAMEFLE PRCELEMS 


In this section several example preblems are presented 
to illustrate the use of the CALICO program in determining 
transfer functions from frequency response data. These 
€xamples illustrate the validity of the kasic concert of 
using the minisization algorithm to determine the transfer 
furction, tut they are by no means an exharstive 
Fresentation of all possible variations of cost functions 
and freguency ranges that say ke employed in determining a 
transfer function representation of a system. For example, 
as will ke shcwn in the example problems, different transfer 
furction representations may be obtained Ey considering 
different ranges of the measured freguency response. That 
is, while the order pf the numerator and denominator are 
fixed by tke user, the parameter values returned frcem the 
prcgram ray vary depending upon the frequency range being 
considered. This may be desirable if seperate low frecuency 
and high frecuency representations are being determined. 
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Also in using the type one cost function, which includes 
both phase and magnitude, the cost surface is different than 
if the tagritude difference alone were used in the cost 
function. This has an effect on the resulting farameter 
values, particularly when lower order transfer functions are 


chcsen tc represent higher order systems. 


The necessary data and the input format employed are 
almcst identical to those used when the prcegram functicns in 
its originally designed compensator design mode. Figure V-1 
illustrates the input data deck used for the first synthesis 
exanple problem. As can be seen the only significant change 
is that the slant has been represented by a unity gain 
transfer function. This has the effect that the resulting 
"ccmpensatcr" represents the transfer functicn modelinc the 


entire syste. 


> 
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Figure V-1 


Transfer Function 


fhe freguency response data, shcwn graphically in 
figure V-Z, was generated from a known transfer function 
using a ccmputer frequency response program. As can he 
seen, twerty magnitude and phase values were selected at 
discrete freguencies over the range from 0.1 radians/second 
to 30 radiansysecond. Examination of the frequency response 
data suggests a transfer function of the fcorm of a corstant 
gain over a g¢clyncmial of some unknown order. This foro 
will yield a constant magnitude at low frequencies anda 
decreasing wagnitude as the frequency increases. Clcser 
exaBinaticn cf the magnitude curve in the range between 2 
and 20 radiansysecond indicates a slope of approximately 60 
Gkyéecade, isplying a polynomial of at least third crder in 
tke denoninatcr. The phase curve extends frem 10 degrees to 
330 degrees over the range cf frequencies ccnsidered. This 
indicates that the denominator may be of fcurth crdeéer or 
higher. As a first guess a third order denominatcr is 
assumed. The results from the CALICO program are shown in 
figure V-2A through V-2H. The initial guess for the 
denominatcr rects and the system gain were chosen as 1, as 
can be seen in figure ¥-2F. Figures V-2A and V-2D show that 
the resulting magnitude curve approximates the measured 
response well at the lower frequencies with a slight 
deviaticn begining tc appear at the higher frequencies. The 
phase plct cf figure V-2B and the plot cf the difference 
between the actual phase and desired phase shcwn in figure 
V-ZE indicate a significant amcunt of deviation frcm the 
measured values in the higher freguency region. This 
Suagests that ancther pole in the transfer functicn is 
necessary tc mcre accurately model the systen. Thus, the 
prcegram ‘as again executed, this time using a fourth order 
dercminatcr. The results illustrated in figures V~2I 
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thrcugh V-2F indicate that this second form of the transfer 
functicn gcre accurately represents the systen. Since the 
criginal transfer function is known in this case, an 
opcrtunity exists to evaluate the performance of the prograr 
in terms cf the accuracy of the transfer function returned 
after minimization of the cost function based on the input 
frequency respense data. This of course would generally not 
be the case when trying to determine the transfer furction 
of an actual system that was to te modeled. fhe parameter 
values returned frcem the CALICO frogram shown in figure V¥-20 
may be ccmpared with the original transfer functicn given 
below. 


126.0 


t(s) = ——_—__—_—___-—_-_-__ - 
(s + 1) (s + 2)(s + 5) (s + 10) 


As can Ee seen, the values of the transfer functicn roots 
returned frem the prcegram for the fourth order denominator 
are in agreement with the values of the original transfer 
function used to generate the frequency resgcnse data. If 
for the system under consideration only the magnitude 
regépense were cf interest, the transfer function ccnsisting 
of a ccnstant gain over a third crder dencminator nicht be 
of sufficient accuracy to model the system. This judgement, 
of course, rests with the user who is familiar with the 
intended use of the transfer function modeling the systen. 


133 


gated ve 


see 


hed 


= f= -- 


a ses == —— 
ae ath ee : 2S SSS ees 


WZ-A oanbty 


T# uny ‘Aouenberq *sa apnytubep 


qP08- 
© 

qP09- 

esuodser tenjoe ,Ssnonutyuoo, --- — 

asuodsez ~Tenjor --- x 

asuodsel paeratsap --- 
APO0Pp- 
qP0c~ 

oT | 

(Ses/spel) 0°OT SB oy . THO 


AP0c 


135 


€Z-A oernbtg 


T# ung ‘Aousenbezq *sa oseud 


> 
} 
o —jo00E- 
rea , 
> 
xX — —.00z- 
2 
asuodsexz [Tenzoe ,snonuTtjuoo, --- — 000T- 
asuodsex jTenjoe --- yx 7 
asuodsert perztsep --- © a 
‘es ee Ty" 0 
(90s/spez) 0°OT O°T eo 
00 


»00T 


136 


© --- desired response 
X --- actual response 


— --- "continuous" actual response 


-300° -~200°. -100% 


a 


Ci’ 


Magnitude vs. Phase, 
Figure V-2C 


ES? 


fy 
) 
‘a 


-10db 


-20db 


=s0d 5 


-40db 


-50db 


-60db 


-70db 


-80db 
Run #1 


Ieee 


° 


100° 


200 


3d 


2d 


“ae . A 1.04 10.0 (rads/sec) 
Saa 
A 
A 

-ldb 

~2db A 

A 
AA 
-3d 


Magnitude Difference vs. Frequency, Run #1 


Figure V-2D 


ise 


T# Ung 


AZ-A 9anbTtg 


‘Kouenberzy 


°SA sdUeTeTJIG eseud 


002 


oO 


009 


u39 


AZ-A eanbty 


T# ung ‘3nd3no0 TeOTTOUNN zoynduoa 


NOTLONNS YDASNVYL YOAVSVIdAID F4L HOS NITLIIIS BHA 
NT aMOTIV 3@° 111M SOUR? BVVIS A1V4 LHOTH ON FACaaN AL Sad AL 
3SVid WAWINIW SHL 39 ST NOTLINNS HZSSMVYL YOLVSNSdWID F4L 
a ee Mi) Oo Gee 0n0"1= 
070 00 3000000°T- 
Ta ta ala re a URN SA SY 0°0 ~~ “00 3000090°I - 
ee LUVd AYVNIDVAI Iuvd Wey ay 
SLOOY HOLVNIVONSO NOLLINAS YSISVVHL HOLVSN3dWII 
me a 90 F000000°T = 0D 2OUCOO*E = -2000790°E = 0_-SNNNADD'T 
oc oe pee a S JD SHIMN4 ONTONSISV NT SIN3ID1SI309 
MOLWNIWONIG NVOTLONNAS HISSNVHL wOLV SNAG) 
i rrr 00 3090099°T 
S 4) _SY3MNd ONIGNIISV NI SLVA19T S309 
ee  SOAVEEWANONOTLONAS HEISNVEL dOLVSN4de99 
00 S000000°T = NIVO NIILINNS YIdSNVEL YOLVSNIDWID 
eee eS ne "99 3090090°T 
RST. Alene | S49 S¥sMod OVvINIISy NI SIN2191 34309 
YOUVNIWINSG NOTLINAS YIISNVYL GSLVSNAdAS ONN 
PRE 6 ta Se - . 00 3090090" 
in : S 40 SH2MId OVIOVIISY NI S1V3I9144209 


YOLVYFIWAN NOTLOINAS 
_ "OO FOCOONOO®T = NIVO NOTLINAS 


6 a ae LdOZ T dWYX4 SISSHINAS NITLI 


USASNVUL JSLVSNS dN DONA 
MSASVV UL GALVSNSdN9DONN 


eereeee weeeetee eg@aeeeernrene fF eererewaee2ee 


NNJ d3dSKVAL --- JULIL 


140 


NZ7-A aanbty 


T# uny ‘3nd3no [Teotzreumn Aaejandwosz 


» 00 5096L99°S = NIVD NOTLINAS AAASNVYUL YIIISNAdNID GA7ZIWTLdO 


oe re ee ee a ee ee eee _- me = tee tees 


2 ee @®FewrFereewewreeereweeeweewe @weaee Weaeeee —“FBF w2eeeewrerwew—wweerwe ee ewrwaese Z@ @ @erZH2 28 = We wsw®Wea- B= — 


0°0 TOS G tol 2— 
ON 3L618S6°T- O07 1289S S > 1— 
eer eg pa Pe ee eat ae OO 32575 S> I ; GO 20965774. {- 


ee x ee ee Luvd AYVNLOVWI bevd ive sdyuv S$LOO" 
YMOLVNIWINAG NOTLONAS ASISNVYL YILISNGdNId) OI7ZIAT1d3 


ee OMIA IGG Ge GOmeC Cra ccam $0 3099009°¢ GO 24b1066°1 


ae ee S 4D SH=NOd ONIGVIOSY NI $iva19144309 
SOLVNIWINSO NOLLONNS WISSNVYL wtagSN3ewOd Ga7Tdlidg 
Se . oS a > “TW Soescricso-2 
ae. ee S 43 S¥y2MU4 ONIUNIDSY NI S1N319144309 
YOLVUSWAN NOLLINNG AZSSNVUL BILOSNAdw2d GIZINTIdD 


ae ee or ee =e 2 = XTdXO9 WOM4 YODD NYUNLTY SIZ FHL 
_19-36906S2°9 = ANIA NJTLINAS LSU ANWINTW SHL 


eS Oe Be 2 e222 Pee @ we es es] Be BF e2eeweerewswewewre ee Vwewree ws Bee Zeewr een ese @ee swe Owe ZBeeeee Be ee = = «= 


ere cee ee ee eee eee SRP oewre reer erewnwre fee @wrewrewreeeree @e ewe wee wwe ewewee Be 8 wwe er wee ene oe 7 2 eee wee @ ee ewe = we wo oe 


1 acdAL JHE StS Gash ad Jl NINN esao saat 


COOOT = sod GAWD SIWIYL JD YAHANN TWVIDL SHI 


ewr@2 22 22 Fe 2s ee ws TZ FMF es ®Feese Fs De Zeew=e © Besser Fes eewrewzaeesww=~H 2 Beas Bea eZwiw Oe Be wee @ ww = @ we we @ @= 


141 


i at ee | 


yal tt ICICI OIC OIAIOUNN 
ewoortToestoce3coo 2909D0C00 


AUIS te hg, WE ct tte 

MPO BCMA os SielelteleKe' 
QOODDSCODRCCOGDSDODDOSCOOSD 
MINES) SES SUC SOoOGee Be wee? 
BLP DVS Qo SA BeiSoVagoeae 
wee SN SNP SULIT LU eC NE PP et NOD 


NOATMMLTSNNALPCPASONTVOAA |: 


oC Re er ee ee ee ee er ee | 


Dt mt IO SU SS et et IOAN AO 


ie er tte he 

; ! i i | : 

t 2 * . 

eed ond Od cond nd nd ot et PIE AIR NI AJAIOIOIOIA] 

QS aOZeasSseqgsscoscac =] 
am, 
COLL tithe Msi a as beats Gocktas Labs tabi bee Abi 
tad TOO NSP ESSN OS SLO DP UV ee QB em 
CDBLCALSAHEAMNMTS™EFDONMEMLAPE AS 
meee Pe md et CN mag Pe EB UN OO Dg St tN DD 
LAO SDMCAIOIORN Det Seci ONS ca 
Nh FOPAILPOaMPIAO MHA MH OOM 
ZOLPAW PLE MOD FMAM FPNNOO 
Teeteteveeteeeeeeeeteee toe e 
Ce et AE PL OD et et eA AISI NIG IOAIOIU™ 
(reeetttd@=~eb at @ get td bt be 


ot . 


cD : 
CADDO we ODO ent et eet ed od Ot et od ot 
mony toc) eo ot COOSlIsT O25 
e2) i | 

Et TOR he Pie Ae [a PO See i ey ee 
ODDP AA ap HOSay GAA DDADS 
OR TIVE SF SAIN SAA MM Met WO 
FE NN A met et ON ee BM CIO CY 
ROT TURE OD ALS TOD AS AOD 
COLON ON LN IN mt VS TOO met PUA OTT et Se pe 
mT OTOP ODTOMS Pate 20 
Mee @ ee Gee eFeeee © © © ee @© eo e 
Lh fom mad met cand ad TN goed at OES DO ee Vie ir Oh 
a) Aes Fe bot poh pon 


{ 
! 


D 
Z 
2 


PROVO COU HHO 0 00 Hatt ete atte 
(ei SALLE > <9 Am Kao Kamel Sees Pla AB J ins Y Sm ae Yet Lae Vm fap | PT a | ee 
= ] 
BBR ONL TS US He LS AR STEED RR IT 
BEPOMO TAM AIO MAN DAMON Qo 
CODED ODO SFO NTANO UD aN ot 
TP ALS GS EDS USN ced ed et mt SO STN SN ee LO ND 
hme YER PAID AIA MH SID AITR LON GS CO me ALE 
meme CEN SY LY med eS od ON) ST PLY PM ate Ue YD 
ZADMBIANTO FROANDYQTNDINEMS 
Dbe3retewneee eee eetrteet @ e & @ 
Let end coed eed ome SONS met OY hm ret et ALTO) PUTT SD 
= at Ne 
| | 

eee oe a = OS DOD [29> 

ee ee 

PION RABE Uy UES LOS ahs UE OU UO OA PRR OR) CLA ICE 
22) DO DSP SI SANS Op Gt'S) Sears) VIST2=) 
ONDADIOWWORIDASTOOODOODO4 
ZRODODOIGCNODR WAS FOTIOOSCOD 
RPO SDT POOoag WMe Swaoe ls] a) 
DONDN OND OT DOSS OAT 
SAPNOST ATI ATT OPK ION YN 
wi © © oe ee Gee ee ee He ee FeO ee Do 
LD met ASS SSRINAL DA POO SRN 
Lo. 


142 


—_— om em em eee ee ee ss ee es ee ee es ee es ee ee ee ee ee ee ee ee ee es oe ee ae oe i oe = oe 


EQUATION INOICATES 


OF TREE CHARACTERISTIC 


Computer Numerical Output, Run #1 


Figure V-2H 


IZ@-A 9anbtgq 


zw# uny ‘Aouenbezqg *sA epnytubew 


Q°T 
(o9s/sper) 0°OT a 
osuodsez Tenjoe ,snonutyjuosd, ---—- — 
osuodsexz Tenjoe --- x 
osuodset peatsep --- © 


qP08s- 


qP09- 


qP0t- 


qP0c- 


aP0 


qpoe 


143 


C7-A eANHTY 


cT# uny ‘Aouonberq *sA oseud 


OCE= 


00c- 


(Ses/spel) 


osuodsez Tenjoe ,snonutyzuos, --—- — 
asuodser yenjoe --- x 
osuodsert peritsep --- © 


144 


© --- desired response 


X --- actual response 
See — “continuous” actual response 


cae)(6 10 =200° -100°3 


Magnitude vs. Phase, Run 


Figure V-2K 


145 


-10db 


=zZ20ab 


=30dr 


-40db 


-50dk 


-60db 


=70db 


=80db 
#2 


100° 


(d90s/spelz) 


T2-A 9aAnbty 


c# uny ‘Aouoenberg *sA soueressIG opnytubew 


qP90°- 


apro°- 


qPc0°- 


APO 


Pz0°* 


qPr0° 


146 


(S0s/spelz) 


WZ-A 9rnbty 


c# uny ‘Aouonbezg *SsA soueTeTJIqG eseud 


(> an) 


147 


NZ-A eanhtg 
cT# Uny ’3Nnd3zNnO TeoOTASUNN Aeynduog 


NOILIWNS BSdSSNVHL YOLVSNGdNI9 34L YOOX NOTLAIIS =HL 
eee MD) eG lec 7 en veneVve Pore OH sags oath fad At 
JSVHd WOWINIWN SHL JO ST NOTLONAY YHASNVaL YOLVSNSdhHUD JHE 
is st, a ees ss 10°C) 00 HOOIVIO°l- t—~*s 
- G0 00 3900990°T- 
a 16. 00 3000099°I- 
RE: Cag ct i fe ht es, : : 0°O ~— Of 300N0D0°T- : 
<i AEE Oe, aT) CR ee a  LdAVd AUYNIOVWI Luvd W3I4 22UV 
S100¥ VULVNTHON SO NOLLONAA vYasSNVed WOLESNIdWOD 
= 00 2000000" I ——_00 100 0000-4 09 3000000°9 00 300090375 2 00 200000021 
7 id Na S 42 SYSMCd ONTIONGISV NI SINII91343909 
YOLVNIWA Na NOI LINNS HAISNVHL DOLVSNIdWII9 co 
a i a a ws Se mi aa atta os sn eS ns a as maaan te SS eo Sg a an amo eae = Je ares treaeen ae ee ea 5 =< 
oo 3000000° I aa 
S 3O SYSMOd INIONIDSY NI SINJID1ISI309 
7 ce POV SWAN NOTLONO A HASSNVYL JOLVSNAGHOD 
00 3000000° T = NIVO VNOILINAS HIASNY UL YOLVSNIdWO9 
- aan —_ - — ~- 00 3000000 I 
ee ee ee S 32 SHY2MOd OVIGDNEISV NI SINS1I14S309 
POLINIWONWZO NOTLINNA YSASVV AL GSLYSNSdNIOINN 
near ee Ss oe See ae 00 5000000°T 
oa a i ne nee te nn S JO SYSMOd ONIONIOSY NI SINIS19144339909 
YOLVYSWAN NITLINASD ASSSNVYL OSLVSNSdNOONN 
EES" FO FOCOONOT I = NIVS NIILINGD HIISVVYUL DAILVSNSdNIDINN . 


ow oe oe om ab Ow ew =e ow om am om om am = 8 wm aw Om ae ow em = aw om aw ae ew em om am 8 OF am SE OP ew ow oe oe oe Oe wm ee Pw ew we ee ie oe 


LdQ2 T dWVX3 SISSHLNAS NOTLONAS AASSNVAL --- 3TVIITL 


O@-A eanbtga 
c# uny ’3ndy3no TeotzASoUWNN AeQyNndUOD 


— 0 JOv909e* T = NIVD NUTLINAS YIASVVYL YILGSVAdNID O3ZINT Id) 


eee — 9g 1)-2311895°6- 
0°0 00 JOLERES*T- 


0°0 00 J0RzyS0°S- 
ee eee 8 O0 JEISTES*S 


AYVNIOVWI Lavd W339 ay 
TLOWNS YsdSNVUL YOLISNAdNId9 O37 INL 
Ble Crea 


df 
que)» TO 3ACSECOT*T OO 301 22c (a0 


TIO, Ra ne, Sa ee UVd 
AILVNIWINAG NS 
wae - £0-396L22T°9 ~=—..00 J6S600T°1 00 46 
SN es own} seme anne an mee mewss-mre Seren S 4) 
YILVNTWINSC NOILOINNS did 


eee Oe we wee ew FeO we Swe ee ee ese eo wee | = 


00 36568TL°2 


—oaw ow ewe ewes Ow eww 8 ee wo ee 2 ee Se SS Swe e Ze Be wwe eZ ow we we we Oe 


ee ES ae Ale OE SE OR = RE a NE I RR Oe Pee me ee ee ee ee ee ee - 


| oe S 40 SHEMC¢d ONTONTOSY NE SINSII 
~ me EVE SWAN NSTLONAS UIISNVYL MOLGSNGdAID D3ZIwTIdd 


$$ ______--—— we em XT KOT ADAA-SGID-NAIVLSY BOUYA FHL 
7 a . - ee ; SO<SEFOTIG°%T = 2AIWA NOILINIS L509 WAWINTW 3H1 


en ere ee Re ee SE a oe OR .- 2 
~~ wee Ow ewe Bw Ow ew ew @e ZF wee 2 SO Se Oe we we Te es Owes Ss ee 


es nn cr emer ee ESM Ce Si) Ta OMe) sala) ea 


Wd QUIWI SIVIdL JD AGHAAN TVLOL IHL 


‘eT ezeren | 


EE ER EE TS DS ES OE SS ES EE TOD SSS ee ee eo ee ee 


149 


tonal t | 
ls ' 
Lt ae md nt meat nd met met od OI NT OG OI IU OUI AIOU AION 
6 - 


CURSO UE 


SBaIeS SOrOoe Pe Soreoeeeye 
[XID DOC MAsorwes =OS9000050 
“ise eu SECON SaASSECE60S05 
(4S Ss 
Rl OM SANT LALIDOANIM MOR SG SF atBO 
MNDMDMNiATeA Data ODD HN 
“1. @ @e @ ee ff @@ &e@e@ © & @ & ee & Oe 
NO DN BN mt et IQ AI QUI VY 

peprptpetoreorrrprrarnrea 
) 


t 
u 
4 . 
' ’ 1 
4 


met. om 


mt end mend md ened gern et carey eet AY AI AI OCA OR AICI NAICS 


ee er) oa eee IO) | 


SOP eeeooesscesoseceoo 


Cay. Ga her ee eae ee Hy i 
ae (Wy WAR SSS SeReuR cosas 

SM mes NOM BOSONS COA PHO 
me UMNO OUI seed ees Ole OO 
: SAW AIT AIA 
DF AO IOIA PUA Pe eee cm PIANO 
LOMA SOD Nat OP HAP ODA AN 
wee a ee er 2 | 
aad ed Be LT et et I AAAS NI OS AS, 

Erp +t eee rea 


—= 
roots] 
TO OO B39 Hh AtO0 OO eh tet Ht et ted 
ee ee ee OS 
(oy i] 
<“biu on Vite, THEI  AIRe, EIR, ali Ded 
S5a5 ow Dein CAM OINAMOD DL 
DAT TADRA SAAD HAaOGAM ANAS DO 
lat NINE SNS Pt mee er Adee CN AIAN 
Sie S) Sipe ate Ss 1 eo a) ON ee San ee) 
CLEVE ORLA md DDD mS FUVLA AO ee Re 
mPOM OTOP TRAST OMNDTAaArEeeroO Oo 
vy? eee @eeeeeteeet*t eet Ge eeee 
ce ad cand ed peg meet TON pred eet OSS ND ed mtd A) A A 
as Cee Ee ke 


ROQDO OO HOOD eww rl td at med ed ont pet 
a, 4 
ES as pers wi Te eee RR 


WIN TONIC NOUR RAMOS eS Sh mas. 


FIAIO ARS NOUS MATUORIOC x 
VOM AN STU MINS AND DOOM SONI 
PDN DS Tet OS DD FP ANIN AG TIM UIE 
PVC NE ES me TT ELEVA DON SU TUS et me pe 
ZPEBOTODVOOS QMO Payer v.g 
Oeeteeeteeteeeet Cee et © @ © @ @ 
ad et ma rag wet Fs ed YN AN AP 
= Peet) et or ea 


‘ ] ; t $ ‘ 
QS att oat ot et IO SONOS OD Dato 
99005593056595995 = ol eo isin ap | 
¢ ¢ 9 bot (4 
Chore seee eo shah Ibigcleks ips sty tere tees 


PLO BOT Oe 
~3S9D073 8 Seog eee eee 
IPOS ISPESIRSE TE SSSSEIS 

SE DBD) 2D) NOP SAME AP Ops FOBVDDSIIO 
SRSs Oar isp Laos POO = DGD ee oe a 
SORA QIOSDIN DESDE ae DAY TIBI 


IRM 1PM PORO SOM YO BDC INO 
or er ee er ee ee ee er 2 a ee Sr 2 ee 
LD mt OD INI N  OS met NX 
th 


150 


ese were fF fe 2 2 ee wPeew ww ewwT ew @®e oF ewes e@wewrewrewerewere Woes ese ewewewereeewrea @ee ws] owe eo |] @ & oo aw o oe & ow ow oe 


ECUATION INDICATES 


CHARACTERISTIC 


nasa ewe eww ee CBD oeews 2 we ST Tew 2 2 | eS! ew | ae ew Sw ow = © ow ew] oe Sw oe oe oe eee oe Ew oF oO oF ow oe ow ee @ o® OF ew ow oe oe ow SO Oa oO oe a oe 


Oo ee + en rere oe ence wee Sr erEnes We cq =eenn come © ry se oes 


wen wwe = - 


Computer Numerical Output, Run #2 


Figure V-2P 


As a second synthesis example, consider the 
freguency response data shown in figure V-3. Here the 
magnitude and phase response for a system are showsn at 
twenty distinct frequency values over the range frem 0.1 tc 
100C.0 radiansyseccnd. As can be seen from this figure the 
positive slcpe of 20 db/decade of the mgnitude curve at the 
lower freguencies suggests a transfer functicn with at least 
one zere. This is also confirmed by the fhase resfecnse, 
which at the lower frequencies initially has a value cf 90 
degrees. At the higher frequencies inspection of the 
magnitude plot shows a negative slope of apffroximately 40 
dktydecade, indicating the dominance of a third order 
dencminatcr in this area. The phase plct approaching -180 
degrees at these higher frequencies alsc suggests this. 
Thus, initially a form for the transfer function consisting 
of a first crder numerator and a third order denominatcr was 
selected tc sodel the system from which the frequency 
response data was obtained. The results are shcwn in 
figures V-3A thrcugh V-3H. As can be seen fron the 
graphical ottput presented in the figures, this particular 
form of transfer function appears to match the _ 2reasured 
freguency response guite well over the entire range of 
frequencies ccnsidered. The rcots and gain cf the transfer 
function rettrned from the program are shown in figures V-3G 
and V-3H. As in the previous example it hapgrens that the 
€xact transfer function used to generate the measured data 
input tc tke CALICC program is known. This transfer 
function is given Fellow in order that the reader may ccmfare 
the values cf the orjginal transfer function parameters with 
these returned frcm the progran. 
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9000s 


t(s) = ————________ 
(s + 2.4) (s + 5) (s + 150) 


As can be seen frem the data presented in figures V~-3G and 
V-2H, the parameter values returned from the program are in 
clese agreement with the original transfer function values. 
While it is urlikely that they will ever be exactly the same 
as the criginal values there is certainly nct a significant 
difference between the two sets of values fcr most 
engineering furroses. 
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fete CHARACTERISTIC EQUATION INDICATES 
Figure V-3H 


Computer Numerical Output 


A frequency response of a system with resonant 
peaks, shcewr in figure V-4, is ccnsidered in this example. 
The same example is considered in ref. 19, where a 
combinaticn cf graphical and analytical techniques are used 
to deterpine the transfer functicn. Initially a transfer 
function ccnsisting of a first order numerator and second 
order dencninator was assumed in trying to mcdel the system. 
This chcice was based on the low and high frequency 
Magnitude curve slopes of approximately +20 db/decade and 
-20 db/decade respectively and the low and high frequency 
phase values cf approximately +90 degrees and -90 decrees. 
This fors failed to match the dip in the magnitude curve in 
the regicn cf 10 radians/second and the jump in the phase 
curve in the same frequency area. After increasing the 
crder of the numerator and denominator, a transfer furction 
consisting cf a third order numeratcr and fourth order 
dencminatcr was found to provide a phase and gain curve that 
closely approximates the measured frequency response. The 
results cf this are shown in figures V-4A through V-4H. The 
nuserical values of the transfer function parareéeters 
returned frc@ the program clesely approximate those of the 
actual system from which the frequency response measurements 
were oktained. The interested reader may verify this fact 
by comparing the values shown in figures V-4G and V-4E with 
these of the actual system given in ref. 19. 
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Figure V-4C 
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In this final example problen a low crder 
approximaticn to a known higher order system is illustrated. 
The original system transfer function under consideration 
consists of a seventh order numerator and an eight order 
dencminatcr given by 


7 6 s ® 3 2 
16s +4E3s +6010s +36380s +122664s 4222088s +185760S+40320 


8 ? 6 s 4 3 2 
Ss +36s +546s +4536s #224495 +67284s 4+118124s +1095845+40320 


An analytical method of approximating this high crder 
transfer function is discussed in ref. 21. Initially 40 
discrete fpcecints, representing the frequency response cf the 
high order transfer function were selected OVEL the 
freguency range from 0.01 radians/second tc 100.0 
Tadians/second. In order to be able tc form a ccmparison 
with the reduced order model obtained in ref. 21, a transfer 
functicn ccnsisting of a first order numerator and a_ second 
crder dercninator was assumed. The results returned from 
the program tsing this form of transfer function are shown 
in figures V-5A through V-SH. It is not surprising that 
there is a significant amount of deviation in the sagnitude 
curve in the vicinity of the resonant effects of the higher 
order system nor that the phase response cf this reduced 
order mcdel dces not match the higher order phase response 
in the hicher frequency range. In order tc compare these 
results with analytical methods for forming reduced order 
ucdels, the frequency response for the reduced order snodel 
of ref. z1 is shown in figures V-5I through V-5M. The 
nugerical values returned from the CALICO prcgram are cf the 
same crder cf magnitude as those found by the analytical 


technigue, rut as can be seen frem the magnitude curves of 
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the two apprcaches (figuresV-5A and V-51), the analytical 
method matches the higher order transfer function extremely 
well at the lower and higher frequencies and makes little 
attempt tc match the higher crder system in the france of 
freguencies at which the rescnant peaks occur. The results 
returned frcemp the program, on the other hand because they 
are a functicn cf the deviaticn of the response frcm the 
desired values over the entire range of frequencies, tend to 
distribute any error between the responses over the entire 
range of frequencies being considered. A slightly different 
representaticn of the system may ke okttained by varyirg the 
frequency rarge under consideration and the type of cost 
function used in the program. For example, figures V-5N 
through V¥-5U illustrate the effects of using a type twe cost 
furcticn and limiting the frequency range from 0.1 to 100.0 
tadians/seccnd. Again the decision as to which acdel is 
best must te made by the user in light of the intended use 


of the recuced order representation. 
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Magnitude vs. Phase, Run #1 
Figure @7-5C 
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Computer Numerical Output, Run #1 
Figure V-5H 
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Magnitude vs. Phase, Analytical Transfer Function 


Figure V-5K 
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Magnitude vs. Phase, Run #2 


Figure V-5P 
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Computer Numerical Output, Run #2 


Figure V-5U 


The investigations carried out in this thesis illustrate 
the feasibility cf automating the classical design of 
compensatcrs by specifying a set of desired values fcr the 
magnitude and fhase response of the open lcop system at a 
nugker cf discrete frequency values over the range of 
interest. The results achieved show that the ccmplex 
gwirimizaticn methcd of M. J. Box is quite carpatle of 
handling the ccnstrained Minimization problen and 
conveniently provides a means of avoiding the numerical 
differentiation cf the cost functicn, required with gradient 
search technigues. The incorporation cf stakility criterion 
in the isplicit constraints of the minimization algcrithn 
alsc insures that the program will not design an unstable 
compensator. Ey specifying the desired respfense in terms of 
the magnitude and phase profiles, ccmmon frequency domain 
specificaticns such as open loop bandwidth, phase margin, 
and gain margin have been incorporated into one ccmmon 
format and tke necessity for specialized routines to check 
these specificaticns have been avoided. Those example 
prcktlems presented for which known results had been 
estaklished Ey other techniques show that the algcrithn 
presented prcvides accuracy that is certainly within the 
necessary tclerances for most engineering applicaticns of 
compensatcr design. 


Also, while not specifically designed for the purpose, 
the basic fremise of the algorithm's Binimizing the 
difference tetween some desired frequency response and the 
frequency respense of the system as calculated by the 
ccmputer prcgram while the transfer functicn parameters are 
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varied has been shcwn to be effective in both determining a 
transfer functicn from measured freguency response data and 
in forming lcw order models of known higher order linear 
systems. In the formulation of lcw order mcdels the choice 
of cost functicn and the frequency range specified may he 
USEC tc emphasize different aspects of the frequency 
response which tke low order model is tc approximate. 


While this werk has shown that the kasic algcrithn 
presented can be used to effectively design series 
ccufensatcrs ktased on frequency response specifications, 
consideratle work remains to be dcne in this area. 
Specifically, an option to select that only real foles and 
zercs be ccrsidered as viable solutions to the compensator 
design algorithm would not be a difficult additicn tc the 
program and should prove quite useful. Also, the 
inccrporaticn of a provision for handling pure time delays 
would not reguire significant modification cf the progran 
ceding and wculd extend the capability of the algcrithm tc 
the soluticn of problems involving transport delays. 
Investigaticn of the method presented should be further 
pursued in terms of multiloof systems and systems invclving 
ccmpensaticn in the feedback loop. It may also prove 
worthwhile tc investigate extending this methcd to nonlinear 
ccmponents using describing functicn technigues. 


In the areas of transfer function synthesis and reduced 
order mcdeling, a more detailed investigaticn of both the 
effects cf different cost functions and considering 
different frequency ranges Should be conducted. In 
addition, it appears feasible to weight the cost furction 
over a farticular frequency range by varying the density of 
the discrete freguencies at which the desired magnitude and 
phase prcfiles are specified. Time did not fpergit a 
detailed investigation of this idea and further 


investigatior intc this area may be aprlicable to the 
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program utilization in both the areas of compensator design 
and transfer function synthesis. 
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APPENDIX A 


DESCRIPTION OF SUBROUTINES 


In this section a brief description and backgrourd of 
various subroutines is presented to aid the interested 
reader in better understanding the internal logic of the 
prcgram. Also these descriptions are intended to facilitate 
the. task of anycne desiring to make mecdificatiors ofr 
additions tc the program in order to handle systems cf a 
configuraticn cther than the unity feedback form assumed 
thrcughout this work. 
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SUBROUTINE CONORM 


The furgese cf this subroutine sub-pfrogram is to 
normalize tte coefficient cf the highest order term to 
unity. Mbis is accomplished in a simple, straight fcrward 
manner by dividing all the coefficients of the polynomial by 
the coefficient cf the highest order tern. Thus a 
Eclynomial given by 


n 
i 
P(s) = a es 
2. 
i1=0 
beccmes cn return from this subrcutine 
n 
an e 
1 Z 
P(s) = ) es 
a 
. n 
1=0 


n 
where a is the coefficient cf the s term. The normalizing 
n 


factor a is also returned to the calling prcgran. 
D 


rt shcvuld be noted that the coefficients of the 
norgalized felynonial are returned under the same variable 
name as tke input coefficients. 


Definiticns cf the input and output parameters are: 


Input Variables 


ee ———— A real one-dimensional array containing the 
ccefficients of the polynomial whose highest crder 
ccefficient is to be normalized to unity. JUpon 
return from the subroutine this array will ccntain 


the ncrmalized coefficients of the felynomial. 


pele 


Ho ~-~-------- The dimension cf the array containinc the 


ccefficients tc be normalized. 


VNOEM ------- The normalizing factor, returned tc the 
calling prcgram, by which all coefficients have keen 
divided. 
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SUBROUTINE PHASE 


This sukroutine sub-program is used to calculate the 
Fhase angle cf a pclynomial at a particular frequency. In 
order tec circumvent the difficulty associated with standard 
arctangent rcutines, which return values limited tc _ the 
range —"Teg-6 2, the phase contributions from the 
individual factors are calculated and then summed to form 
the total phase angle at a particular frequency. That is, 
if the pclyncnial is factored and represented as the product 
eof first crder terzs of the form 


n 


P(s) = ial (s +z.) 
1 


i=0 


where in general z may be a complex number of the form 
i 
Zz =a._+dgjk, then at a particular frequency s = je the 
i af i 


phase contrikution of the individual factors may be 


evaluated as 


(bb, tw ) 
Z 


a e 
a 
and the tctal phase angle at the particular frequency is the 
¢um of these individual phase contributions. fhe phase 
contributicn cf each of these first crder terms will be 
within the principal values of the arctangent routine cf the 
ccaputer and the problems generally associated with 
determinaticn cf the phase for higher than first order terms 
are avoided. 


Definiticns of the input and output parameters are; 


Input Variables 


, were crr--- A reak one-dimensional array containins 


real sarts cf the rocts of the polyncspial. 


ee A real one-dimensional array containing 


isaginary parts of the roots of the felynomial. 


cietentntetentetentetmeent An integer value specifying the number 


rects cf the polynomial. 
EGA -----~- The frequency in radians/sec at whick 
value cf the phase is to be computed. 
Output Variables 


coo ----——— The total phase contrikution of 


pclyncnial civen in degrees. 
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the 


the 


of 


the 


the 


SUBROUTINE PINVRET 


The fpurpese cf this sukroutine suk-fregram is to 
rearrange tke coefficients of a polynomial in inverse order 
from that fed into the subroutine. Consider for example a 
Felynomial arranged in descending powers of s 

E(s) = a a + a ae + c200) todays 1 a 
D n-1 a 0 
Upen return frcm this subroutine the coefficients wculd be 
arranged in ascending powers cf s. That is the pclyromial 
wculd be arranged as fcllows 
E(s) =a +¢a4S + «ee + 48 Pe + a S 
Ad i n-1 n 
While kcth representations are mathematically equivalent 
the sequence cf the coefficients must be taken into acccunt 
when using various pclynomial manipulation techniques tsithin 
the computer. The method used in reversing the sequence of 
the input ccefficients consists of temporarily storing the 
first coefficient and sequentially shifting the remaining 
coefficients. The coefficient in temporary storage is then 
Flaced in the last coefficient location. This process is 
continued with the use of nested do loops and each time the 
shifting operation occurs the coefficient in temporary 
storage is laced one location lower than the previcus 
number retrieved from temporary storage. 
\ 

This stkroutine requires no cther subroutines or 
functions tc carray out the procedure. The size of the one 
dizensicnal array of coefficients that may te rearranged is 
determined ky the dimension cf the array in the calling 
progran. 


Definiticns of the input and output parameters are: 
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Input Variables 


ae -------- A real one dimensional array containing the 
ccefficients tc be rearranged. Upon return frer the 
subroutine this array contains the input coefficients 


in reverse seguence from that entered. 


_----------- The dimensicn of the array containing the 


ccefficients tc be rearranged. 
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SUBROUTINE PMPY 


The purpose of this subroutine suk-progragzp is to 
multiply the ccefficients of one polynomial X(s) times the 
coefficierts cf anecther polynomial Y(s) to obtain the 
coefficients cf the resulting polynomial Z(s). That is, if 
it is desired to multiply two polynomials of the forn 


& n-1 n 
a *+astas oe) teat s + aé 
re) 1 2 


P (s) 


and 


© al n 
bo +¢+bst+ bs + «ee + D s + | mn 
0 a € n-1 n 


P(s) 


tcgether in crder to obtain a third resultant polynomial of 
the general fcrm 
2 i-1 1 


P (s) =c t#+6S8S+4+CS *#¢ coe tC 
3 ) 1 2 


if) 
+ 
‘@) 
tn 


i-1 i 

this stkrottine is used in acccmplishing this. The order 
of the resulting pclynomial is calculated by forming the sum 
of the orcers cf the two input pclynomials as 

l1=n +o 

The varicus ccefficients of the resulting pclynomial, which 
are returned ir ascending powers of the independent 
variable, are calculated as the sums of products of the 
input ccefficients by means of two nested do-loops. The 
arrangement cf the two do-loors is such that the approfriate 
product terms are summed in forming the respective 
coefficierts cf the resultant pclynomial. All input and 
output ccefficients must be real numbers. The size cf the 
Eclynomials that say be multiplied are contrcelled by their 
respective dimension statements within the calling prcgram. 


Definiticns of the input and output parameters are: 
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Input Variables 


X 22-2 --- = A real one dimensional array containing the 


ccefficients of the first polynomial tc be 
multiplied. These must be arranged in ascending 
powers of the independent variable. 


------- The dimension cf the array X c ontaining the 


ccefficients of the first polynomial. 


-—----- A real one dimensional array containing the 


ccefficients of the second polynomial tc be 
multiplied. These must be arranged in ascending 
pewers of the independent variable. 


ee is The dimension of the array Y containing the 


ccefficients of the seccnd polynomial. 


Output Variables 


==e---- A real one dimensional array containirg the 


ccefficients cf the resultant pclynomial. These 
ccefficients are arranged in ascending powers cf the 


independent variable. 


------- The calculated dimension cf the array Z 


ccntaining the coefficients of the resultant 
pclyncpial. 
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SUBROUTINE PVAL 


The purgpese cf this sukroutine sub-frcegram is to 
evaluate a pelyncnzial having the complex variable s = o + jw 


and real coefficients a. For a given polyncmial of order n 
2 


of the fcrao 


i=0 


s is declared as a complex variable within the subrcutine 
and the indicated operations of summation and multiplication 
are carried cut by a simple lcop operation taking advantage 
cf the asscciative and distributive laws cf multiplication 
and additicn. This subroutine requires nc other 


suk-programs tc perform its function. 


Definiticns of the input and output parameters are: 


Input Variables 


2 ————--—-—— A rveal one-dimensional array containing the 
ccefficients cf the polynomial that is tc be 
evaluated. The coefficients are assumed to be 
arranged in ascending order according to the fowers 
of s with zero explicitly input for the coefficient 
of any missing powers of s occuring in the 
pclyncazial. The maximum size of this array 1s 
ccntrclled by the dimension of the array withir the 
main frogram. The values of the coefficients in this 
array remain unchanged after execution cf the 
sukroutine. 
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= —————=-- The order of the polynomial tc be evaluated. 
This value is one less than the dimensicn of A(I). 


a ————— The value of the real fart cf the ccnuplex 
variatle s at which the polyncnmial is to be 
evaluated. 

a ————— = The value of the imaginary part cf the 
ccuplex variable s at which the polynomial is to be 
evaluated. 


Output Variables 


a ————— -- The real part of the complex resultant value 
of the polynomial evaluated at s = o + jw, 


————————~——= The imaginary part of the ccmplex resultant 
value cf the polynomial evaluated at s = ¢ + jo. 
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SUBROOTINE ROUTH 


fhe furgcese cf this subroutine sub-frogram is to 
determine if a given input pclynomial possesses right half 
Flane rocts ty using the Routh stability critericn. The 
EClynomial ccefficients are input in descending pewers of 
the independent variable. From these coefficients the Routh 
array is calculated and the first column of this array is 
searched for an indication of the existence of fright haif 
Flane rocts. Thus from an input polynomial cf the forn 


n n-1 
E(s) = bs +b Ss + cee + bDStd 
n Do 1 ) 
tke Reuth array 
st 
1 row b b L ete os 
n n-2 n-4 
nd 
Low b b b siSraveeve 
n-1 n-3 n-5 
rd 
Ss Low a a a bre ana esa 
31 32 33 
th 
4 row a a a alereteieiers 
41 42 43 
th ° ; 
n row a a a arerevetee te 
n1 n2 n3 


st 
(n+ 1) rCcw a a a 
(n+1)1 (n+1) 2 (n+1) 3 
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is formed, where in general the entries from the third row 


cnward are given by 
[2 (n> 1} il? (n-2) (i+1)] 7 [2 (n-2) il (n-1) (i+) 


a 
(n-1) 1 
The first cclumn of this array is then tested to determine 


Cs 


if any of tke elesents are less than zero. Note that since 
all coefficients are limited to positive values by the 
search area restrictions in cther parts of the prcegram no 
prcvision has been made to handle the case cf the 
coefficients cf the polynomial being all negative. Thus a 
pclynomial with all negative coefficients should be 
aultiplied ky -1.0 prior to using this particular 
subroutine fer testing for right half flane roots. 


Definiticns of input and output parameters are: 


Input Variables 


Y 2-2-3 ------ A real one dimensional array containing the 
ccefficients of the polynomial in descending fewers 
of the independent variable. 


N ----------- The dimension of the array Y. 


Output Variables 


a3 >-———— Integer output indicating the absence or 
presence of right half plane roots 
O-~-Indicates the polynomial has no right half 
Flane rpots. 
1---Indicates the polynomial has fright half 
Flane roots. 
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SUBRCUTINE SEMEL 


The furpese cf this subroutine sub-frogram is to 
calculate the ccefficients of a polynomial from the rects of 
the polyncaial. The calculation is done using ccmplex 
arithmetic in order to te able to compute the polyncmial 
coefficients when complex roots are present. The resulting 
EFclyncmial ccefficients are, however assumed to he real. 
Thus, if ccmplex roots are present they must be in conjugate 
fairs. Alec the coefficient of the highest order term is 
set egual to unity. Consider a polynomial of degree n 
having n rocts. This gay be expressed as 


n 
P(s) = | | (Ss + 2.) 
1 

1=0 


where in general the z's may be complex. The coefficients 
1 


of the polyncmial that results when the above indicated 
multiplicaticn is performed may be computed Ly the fcruaula 


n71 
A. = (~1) e > (product of roots taken n-i at a time) 
an 


Thus the coefficients of the pclynomial are given by 
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a = Z @Z @Z @®Z #®...°2Z 
1 2 3 4 


n 
a = 2 @©2 @#Z #22492 + Z @Z @#Z @©...°2Z 
1 12 3 net 23 4 n 
a = Zz + Z +7 Z + Z 38 4 Tr lsie 6) Te 
a = 1 
D 


The largest crdeyr polynomial that this particular sukrcutine 
th 
will handle is a 20 order. If larger pclynomials are to 
be considereé the dimension cf R within the subroutine must 
be changeca. 
Definiticns of the input and output parameters are: 
Input Variables 
a --—-———= A real one-dimensional array containing the 
real farts of the polynecnmial rocts. 


= ————— A real one-dimensional array containing the 


imaginary farts of the polynomial rccts. 


a — An integer specifying the order of the 
pclyncpial or equivalently the dimension of the 
arrays containing the real and imaginary parts cf the 
ECCtS,. 


Output Variables 


a — A real one-dimensional array containing the 


resulting ccefficients cf the polyncmrial arranged in 
ascending powers of the independent variatle. The 
dinension of this array must be one larger than the 
digensicns of the arrays containing the real and 
imaginary parts of the roots. 
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APPENDIX B 


PROGRAM LISTING 
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APPENDIX C 


CESCRIPTION OF FORTRAN VARIABLE NAMES 
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A (I) 


ADAMG (I) 


AMAG (I) 


ANGDEL 


ANGDER 


ANGNUL 


ANGNUE 


E (L) 


CO (I) 


COF (I) 


COFASN (I) 


UNIT 


None 


db 


db 


Legrees 


Radians 


Tegrees 


Radians 


None 


None 


None 


None 


DESCRIPTION 


Arra ood the coefficients 
of he plant ransfer function 
Numerator . 


Array containing the desired crfen 
loop magnitude respense in db at 
the discrete frequencies specified 
by the designer. 


Output array Containing the 
magnitude of the cren Loop 
frequency response in db at. tte 
discrete frequencies specified. 


Phase angle contribution of the. 
open loop system transfer function 
denominator at a specific 
frequency. 


Phase angle contribution of tke. 
open Poae system transfer function 
denominator at a specific 
frequency. 


Phase angle contribution of the. 
open loop system transfer function 
numerator at a specific frequency. 


Phase angle ccntribution of the. 
open loop system transfer function 
numerator at a specific frequency. 


re! containing the coefficients 
of the plant transfer function 
denominator. 


One dimensional working array used 
primarily in reading POLI REG onl 
coefficients and infuting data to 
root finding subroutines. 


One dimensional werking array used 
primarily to store resulting 
coefficients used in checking 
accuracy of root finding 
subroutines. 


ee containing the coefficients 
of the cpen loop transfer function 
numerator. 
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FORT 
VARIA 


fest>s] 
te 
(xl oe 


COFASL (I) 


COFWRK (I) 


DELPHI 


DELPSI 


DIFF1 (I) 


DIFF2 (1) 


DMAG (T) 


DCPHSC (I) 


DWFQ 


FREQ (I) 


GAIN 


GAINCC 


UNIT 


None 


None 


Degrees 


Degrees 


db 


Tegrees 


None 


Degrees 


Radians 


None 


None 


None 


DESCRIPTION 


Array containing the coefficients 
of the cpen loop system transfer 
function denominatcr. 


One dimensional werking array used 
primarily to temporarily stcre 
ag Tae coefficients in 
performing polynomial 
manipulations. 


Difference between two numerator 
phase angle ccntrikutions at. 
Successive discrete frequencies. 


Difference between two denominator 
phase angle contrikutions at, 
SuccesSive discrete frequencies. 


Difference between the desired and 
actual magnitude of the open 2088 
frequency response at the discrete 
fe auency, values specified in the 
input. 


Difference between the desired and 
actual phase of the open lccp 
frequency response at the discrete 
fpegzency values specified in the 
input. 


Array containing the desired 
magnitude response at the discrete 
frequencies specified by the 
designer. 


Input array containing the desired 
phase profile cf the open loor . 
System at the discrete frequefcies 
chosen. 


Interval to be used if discrete 
sec ae cy values are to be 
incremented linearly. 


Logarithm of the discrete : 
frequencies used in constructing 
abscissa for Bode frlots. 


Gain of plant transfer functicn. 


Gain of compensatcr transfer 
function. 


NS). 


IDB 


IERBXE 


INFORB 


IPLOT 


IZOH 


KNOW 


NEODE 


NCD 


NCN 


NICHCL 


NMINFS 


NOMEG 


NSD 


UNIT 


None 


None 


None 


None 


None 


None 


None 


None 


None 


None 


None 


None 


None 


DESCRIPTION 


Input integer to be set egual to 
unity if the desired open lcoor 
Magnitude response is input in db. 


Output integer indicating the . 
completion code of the mininization 
routine. 


Alphabetic input indicating if 
coefficients are te ke read ir 
factored cr polyncmial forn. 


Integer input, to ke set to unity 
1£ CALCOMP plots of results are 
desired. 


Integer input set egual to eno! if 
there 1s a zero order hold in the 
error Sa path prior to the 
plant and compensatcr. 


Integer input indicating if 
discrete frequency values are to be 
read in or incremented linearly 
from the minimum freguency. 


Integer input to re set equal to 
unity 1f£f Bode plot is not desired. 


Order of compensatcr transfer 
function denominatcr. 


Order of compensatcr transfer 
function numerator. 


Integer input to ke set equal tc 
unity if Nichol's flot is not 
desired. 


Integer input to ke set equal to 
See es a nopgaininum phase 
solution is to be allowed. 


Integer Sore Eee the total number 

of discrete frequency points teing 

considered in the input frequency 

be 800 profile. This numEer must 
e . 


Order fo the system open lccp 
transfer function denominatcr. 
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NTIFN 


NYQST 


N3 


N4 


CDENT 


CDENR 


CMAG (TI) 


CNUMI 


CNUMR 


CPHSD (1) 


UNIT 


None 


None 


None 


None 


None 


None 


None 


None 


Nene 


None 


None 


LTegrees 


DESCRIPTION 


Order of the system open lccp 
transfer function numerator. 


Order of plant transfer functicn 
denominator. 


Order of plant transfer functicn 
numerator. 


Integer input to ke set equal tc 
unity if Nyquist pict is not 
desired. 


Dimensicn of array containing plant 
transfer function numerator 
coefficients. 


Dimensicn of array containginc 
plant transfer function denominator 
coefficients. 


Value of the immaginary part cf the 
systen epee Zoe Canster functicn 
denominator evaiuated at a 
particular frequency. 


Value of the real fart of the |. 
systen open Sook transfer function 
denominator evaluated at a 
particular frequency. 


Output array containing the 
magnitude of the frequency resfense 
computed for the oren ree te ste en 
transfer function at the discrete 
frequencies specified. 


Value of the pees lnaLy part cf the 
System open loop transfer function 
numerator evaluated at a particular 
frequency. 


Value of the real part of the | 
System open loop transfer function 
Numerator evaluated at a particular 
frequency. 


Output array containing the phase 
of the frequency response of the 
open loop System at the discrete 
frequencies specified. 


Zot 


CPHSR (1) 


SAMAG (I) 


SAPHST (I) 


WFREq¢1) 


WMAX 


WHIN 


WS 


XS (T) 


UNIT 


Radians 


db 


Degrees 


Secs 


Radians 


Radians 


Radians 


Radians 


None 


DESCRIFTION 


Output array containing the phase 
of the freguency respfonse of the 
open loop system transfer function 
at the discrete frequencies 
specified. 


Array containing the magnitude of 
the Simulation of the ES eal: 
Optimized system transfer function. 


AruAy containing the phase cf the 
Simulation of the a . 
Optimized system transfer function. 


se naeag period if a zero crder 
hold is present in the systen. 


Array ccntaining the discrete 
freguencies specified by the 
designer. 


Input = peered the maxinun 
frequency value of the range cf 
frequencies being ccnsideréd. 


Input specifying the mininun 
frequency value of the range cf 
frequencies being ccnsidered. 


Output giving the sampling 
frequency if a zero crder hold is 
present in the system. 


Array in which values to be varied 
by the minimizaticn routine are 
stored. At the start this Se 
contains the initial ques of the 
coefficient values and at . 
completion the areas contains the 
coefficient values that yield the 
Minimum cost functicn. 
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